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ABSTRACT 


Mechanical  properties  of  advanced  composite  lamina  are  identified  for  better 
mathematical  modeling  of  composite  laminate  for  structural  analysis.  Each  lamina  is 
treated  as  an  orthotropic  material  under  plane  stress  state  and  are  assumed  to  be 
transversely  isotropic.  Four  stiffness  properties,  (Eh  E2,  n12,  Gi2),  are  treated  as  design 
variables  for  minimization  of  a  performance  index.  The  differences  between  analytically 
obtained  and  experimental  natural  frequencies  for  the  specimen,  along  with  a  proper 
weighting  scheme  for  each  mode,  are  minimized  using  the  optimization  routine, 
‘fmincon’  in  the  MATLAB®  optimization  toolbox.  The  modal  assurance  criterion  is 
utilized  to  construct  the  weighting  to  express  the  degree  of  correlation  between  mode 
shape  vectors  obtained  experimentally  and  derived  analytically. 

This  study  requires  a  series  of  experimental  results;  natural  frequencies  and 
corresponding  mode  shapes  of  the  specimen.  A  computational  tool  has  been  developed  as 
a  result  of  this  study.  Numerical  examples  are  investigated  to  demonstrate  the 
performance  of  this  approach.  Further  study  with  experiments  may  show  practical  benefit 
of  current  method  for  characterization  of  mechanical  properties  of  advanced  composite 
materials. 
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I.  INTRODUCTION 


Based  on  the  profound  benefit  of  composite  materials  such  as  high  specific 
stiffness/strength  and  toughness,  etc.,  it  has  become  common  to  utilize  composite 
structures  in  many  applications.  Laminated  construction  is  the  most  popular  type  of 
composite  structures.  Owing  to  many  researchers’  work,  the  advanced  composite  material 
can  be  mathematically  modeled  under  well  defined  assumptions  [1].  A  lamina  is  the  basic 
building  block  of  laminated  composite  structures.  Five  independent  constants  in  stress- 
strain  relationship  for  transversely  isotropic  materials  are  reduced  to  four  when  we 
assume  the  plane  stress  states  [2]. 

As  material  constants,  the  four  independent  mechanical  properties  should  be 
determined  experimentally  to  be  used  to  model  the  composite  structures  mathematically. 
There  has  been  extensive  research  to  characterize  these  lamina  properties.  A  series  of 
static  tests  for  several  specimens  have  been  proposed  since  1960’s  [3].  A  number  of 
researchers  have  tried  to  get  these  properties  from  a  dynamic  test  of  a  specimen  [4-10].  It 
is  based  on  the  fact  that  the  measured  modal  properties,  e.g.,  natural  frequencies  and 
mode  shapes,  are  functions  of  physical  properties  of  the  structure.  Inverse  problems  have 
been  solved  to  find  the  parameters  in  the  mathematical  model  which  can  match  the 
analytical  modal  properties  with  those  of  the  real  structures.  As  an  analysis  tool,  various 
methodologies  have  been  applied,  including  for  example,  the  Rayleigh-Ritz  technique 
[4~6]  and  the  finite  element  method  [7]. 

Most  studies  have  focused  only  on  the  magnitude  of  the  natural  frequencies  of  the 
analytical  model.  The  smallest  (first)  natural  frequency  of  analytical  model  is  compared 
with  the  first  frequency  from  experiment,  and  so  on.  Natural  frequencies  alone,  however, 
cannot  represent  satisfactorily  the  dynamic  behavior  of  a  system.  Suppose  a  set  of 
specific  mechanical  constants  yield  the  natural  frequencies  which  coincide  with  the 
corresponding  experimental  frequencies.  We  should  confirm  whether  the  analytical  and 
experimental  mode  shapes  match  in  addition.  If  the  corresponding  mode  shapes  of  the 
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analytical  model  are  similar  to  the  experimental  ones,  the  mechanical  constants  involved 
can  be  said  to  describe  the  system  properly. 

In  this  study,  the  four  mechanical  stiffness  constants  of  a  composite  lamina  are 
estimated  through  minimization  of  the  perfonnance  index,  which  includes  the  similarity 
between  the  experimental  mode  shapes  and  the  analytical  mode  shapes  as  well  as  the 
differences  in  the  corresponding  natural  frequencies.  Differences  between  the  natural 
frequencies  from  experiment  and  those  from  analytical  model  for  corresponding  modes 
are  weighted  by  factors  based  on  the  concept  of  modal  assurance  criterion  [11].  The 
weightings  for  each  mode  express  the  degree  of  correlation  between  the  experimental 
mode  shapes  and  the  analytical  ones. 

The  analytical  model  of  the  specimen  makes  use  of  the  classical  laminate  theory 
(CLT)  and  Reissner-Mindlin  plate  theory.  The  finite  element  method  employs  the 
isoparametric  nine-node  plate  element.  The  specimen  is  a  laminated  plate  with  an 
arbitrary  but  known  lay-up  in  a  cantilever  plate  configuration.  A  proper  number  of  finite 
elements  is  used  to  model  the  specimen.  Eigenvalues  and  eigenvectors  up  to  the  fifth 
mode  are  calculated  to  compare  with  the  experimental  results.  The  vertical  displacements 
at  the  center  node  of  each  element  comprise  the  mode  shape  vector  for  comparison  with 
the  experimentally  obtained  mode  shape. 

Performance  index  minimization  is  performed  using  the  optimization  routine, 
fmincon.m’  in  the  MATLAB®  optimization  toolbox.  Four  elastic  constants 
( El,E2,vn,Gn)  have  been  treated  as  design  variables  for  minimization.  During  the 

minimization  process,  the  four  design  variables  are  updated  such  that  the  resulting 
analytical  responses,  i.e.,  natural  frequencies  and  mode  shapes,  match  to  the 
corresponding  experimental  ones.  The  sensitivity  of  the  natural  frequencies  with  respect 
to  the  design  variables  is  investigated. 

This  study  requires  a  series  of  experimental  results;  natural  frequencies  and 
corresponding  mode  shapes  of  the  specimen.  A  computational  tool  has  been  developed  as 
a  result  of  this  study.  All  the  procedures  are  coded  using  MATLAB®.  Numerical 
examples  are  investigated  to  demonstrate  the  performance  of  this  approach.  Further  study 
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with  experimental  results  may  show  practical  benefit  of  the  current  method  for  the 
characterization  of  mechanical  properties  of  advanced  composite  materials. 


3 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


4 


II.  SYSTEM  MODELING 


A.  EQUATION  OF  MOTION  AND  FINITE  ELEMENT  FORMULATION 

The  equation  of  motion 

V*cr  +  /  =  /?v  in  V  (1) 

can  be  expressed  as  Eq.  (2)  by  use  of  principle  of  virtual  work  and  the  divergence 
theorem. 

J  p(5v)T'idV  +  J  (5s)r  crdV  =  J  (, SvffdV  +  J  ( SvfpdS  (2) 

Here  we  assumed  the  strain-displacement  relation  as  Eq.  (3)  for  small  strain. 

£  =  I(Vv  +  Vvr)  (3) 

The  stress-strain  relation,  Eq.  (4)  is  applied  to  Eq.  (2)  to  obtain  the  weak  form  of  equation 
of  motion,  Eq.  (5),  where  we  deal  with  the  case  of  no  external  forces. 

c t  =  Ds  (4) 

J  p(Sv)T'idV  +  J  (Ss)T  DsdV  =  0  (5) 

Introducing  the  displacement  interpolation  ( N ),  the  strain  interpolation  function 
(. B ),  and  the  nodal  displacement  (U),  we  obtain  the  discretized  finite  element  equation  of 
motion,  Eq.  (6). 

MU  +  KU  =  0  (6) 

where,  M  is  mass  matrix  as  Eq.  (7)  and  K  is  stiffness  matrix  as  Eq.  (8). 

Ne 

M  =  YJPNTNdv  (V) 

Ne 

K  =  YjBTdBcIV  (8) 
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B.  NINE-NODE  REI S SNER-MINDLIN  PLATE  ELEMENT  (C°) 


The  laminated  composite  plate  is  modeled  using  the  nine-node  Reissner-Mindlin 
plate  elements  based  on  the  classical  lamination  theory  [2]  and  the  first  order  shear 
deformation  theory  [12].  Each  node  of  the  finite  element  has  three  degrees  freedom, 
( </)x,(/)2,w ),  which  corresponds  to  two  rotation  angles  and  one  vertical  displacement. 

Figure  1  shows  the  relationship  between  the  shear  strain  and  the  derivative  of  vertical 
displacement  in  the  xz-plane  for  the  first  order  shear  defonnation  theory.  Transverse 
shear  strains  are  expressed  as  follow: 


r„ : 


dw  , 


dw  , 

'y — & 
dy 


(9a) 

(9b) 


Linear  displacements  along  x-,  y-,  and  z-direction,  can  be  expressed  as 


u(x,y,z)  =  -z</>](x,y) 

(10a) 

v(x,y,z)  =  -z<f>2(x,y) 

(10b) 
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w(x,y,z)  =  w(x,y) 


(10c) 


And^ ,  (j)2 ,  and  w  are  interpolated  as  Eqs.  (11), 
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f*i 

1=1 

(11a) 

% 

-sT 

ll 

(lib) 

-sT 

Ji 

ll 

St 

(He) 

where,  hj  are  isoparametric  interpolation  shape  function  (Eqs.  12)  for  the  element  shown 
in  Fig.  2  [13]. 


1 

S  =  1 

l  5 

% 

1 

i 

4 

8 

*7  3* 

9  6 

1  4 

P 

r 

k  — - 

- 1 

f - 1 

9  1 

V - 

r=-l 

r=  1 

1 

5  2 

i 

i - < 

t - 4 

i 

s=-l 

Fig.  2  Nine-node  Isoparametric  Element 


hl(r,s)  =  ^(r2 -r)(s2 -s) 

(12a) 

h2(r,s )  =  -^-(r2  +  r)(s2  ~s) 

(12b) 

h3(r,s )  =  -^-(r2  +r)(s 2  +-s) 

(12c) 

h4(r,s )  =  -^-(r2  -fX-s2  +  •?) 

(12d) 
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h5(r,s)  =  ^(\-r2)(s2  -s) 

(12e) 

h6(r’s)  =  ^(r2  +r)(l~s2) 

(12f) 

h7(r,s)  =  ^(\-r2)(s2 +s) 

(12g) 

K(r,s)  =  ^(r2  -r)(l-s2) 

(12h) 

h9(r,s)  =  (\-r2)(l-s2) 

(12i) 

Displacement  interpolation  within  a  finite  element  is  expressed  in  matrix  form  as  Eq.  (13). 


u 

1 

0 

0 

V 

>  = 

...  o 

-zfy 

0  ••• 

w 

0 

0 

-h, 

i  =  1  ~  9 


(13) 


The  strains  components  of  interest  are  expressed  in  terms  of  displacement 
interpolation: 


— =  -zM  =  _zy®^<» 

dx  dx  ~~t  ox 


e  >=-z8t  =  -.y%> 
"  8y  dy  t 


=  + = -Z(M+M) = 3*® 

3v  dx  5  v  dx  [  dv  dx 


dw  x  7  /  (/)  S/z 

rx ;  =  - —  <t>x=-YjhA  +  Z  ^ 

uX  ;-=j  C/X 


5w  ^  x-1 ,  ,  (i)  dht 

/vz  =  T - ^2  =  “Z  hA  +  Z  ^  W/ 

dv  /=1  /=i  dy 


(14a) 

(14b) 

(14c) 

(14d) 

(14e) 
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Strain  components  are  grouped  into  in-plane  and  transverse  components  as  follows: 


Y: 


xy 


rx: 


Y yz 


dh, 

0 

0 

dx 

••  0 

dh, 

o...  , 

A10 
<t>2 2(i) 

dh, 

dy 

dh 

0 

w(i) 

dy 

dx 

-h 

0 

dh 

4(° 

i 

0 

~h, 

dx 

dht 

•  A® 

dy 

W 

zB,{ue], 


i  =  1~9 


i  =  1~9 


(15a) 


(15b) 


The  element  stiffness  matrix,  Ke  is  the  sum  of  the  bending  stiffness  matrix,  if/ 
and  the  transverse  shear  stiffness  matrix,  Kse .  The  shear  correction  factor  is  assumed  as 
k  =  5/6 . 


Ke  =  K/  +  k  ■  Kse 


(16) 


Detail  expressions  to  calculate  the  bending  stiffness  and  transverse  stiffness  are  shown  in 
Eq.  (17)  and  (18)  below: 


K/  =  JJJ  {-zBI  )T  Q,  (-zBj  )dV  =  JJ  B, 


=JJ 


aA  o 

dx 


3 \ 

dy 


0  VK 

dy  dx 

0  0  0 


z 

J  QjZ2dz 


Dn 

D\  2 

Dl6 

Dn 

D22 

^26 

,A« 

^26 

D66 

BjdA  =  JJ  Bl'DBBldA 


dx 


0  — t- 


dh, 
dy 
dh  dh, 


0  0 
0 


i  i 


dy  dx 


0 


(17) 


dA 
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t 


K/  =  j]J  B/QsBsdV  =  JJ B /  }  Qsdz  BsdA  =  J| B/DsBsdA 

t 

~2 


-h,  0 

0 

d \  dfy 

dx  dry 

where  the  off-axis  stiffness  of  lamina,  Q’s  and  laminate  stiffness  matrix,  DB  and  Ds  will 
be  described  in  detail  in  later  section. 

As  the  shape  functions  are  expressed  in  terms  of  natural  coordinates,  we  need  the 
relationship  between  the  physical  coordinates  and  the  natural  coordinates  to  obtain  the 
derivatives  of  shape  function  with  respect  to  the  physical  coordinates.  Using  the  chain 
rule  in  Eq.  (19), 

d  dx  8y  d 

dr  _  dr  dr  dx 

d  dx  dy  d 

ds )  L  ds  ds  J  [  dy 

The  derivatives  of  the  shape  functions  with  respect  to  the  physical  coordinates  are 
expressed  as  Eq.  (20). 

dhj_  f  dht 

dx  ,i  dr 

dhi  dhi 

dy)  l  c!s 

The  element  mass  matrix  is  expressed  in  a  similar  way  as  in  Eq.  (21).  Mass 
throughout  the  element  is  assumed  constant  and  the  mass  density  is  p  . 
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Me  =  JJJ  pNTNdV 

-zht  0  0 

=  \\\p  0  -zht  0  ••• 

0  0  k 


h  0  0  f  —  f — z  1  If  k  0  0 

=  JJ p  0  ht  0  «  J  — z  [-z  -z  l] dz  >  •  ■  •  0  ht  0  ■■■  dA 

0  0  h,  [4  1  J  J  [  0  0  ht 
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=  —\\  0  ht  0  0  1  0  •••0  h,  0  •••  dA 

12  0  0  h,  o  o  12  [  0  0  ht 

i  JL  (21) 

Gauss-Legendre  quadrature  is  utilized  to  integrate  the  polynomials  for  the  finite 
element  matrix.  Three  points  integration  gives  the  exact  value  for  the  shape  function  used 
in  this  study.  For  the  element  mass  matrix  and  the  bending  stiffness  matrix,  three  by  three 
point  integration  is  performed.  To  avoid  shear  locking,  two  by  two  point  integration  is 
chosen  for  the  transverse  shear  stiffness  matrix  [14]. 

Eqs.  (22)  give  the  expression  for  the  matrices  of  the  finite  element. 


=  WW^J) 

i=\  j=\ 

(22a) 

Ks-  =  t Z j B/DsBs}  W,  W] det(J) 

'=1  7=1 

(22b) 

m'  =  pYL  k  w> del(/) 

i= 1  2=1  ',S‘ 

(22c) 

-zh,  0  0 

0  -zhj  0  •••  dV 

0  0  ht 
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where,  r,  and  .v,  are  Gauss-Legendre  integration  points  in  the  r-  and  5-direction, 
respectively;  W,-  and  Wj  are  the  corresponding  weightings. 


C.  LAMINATE  STIFFNESS 


For  transversely  isotropic  lamina,  there  are  five  independent  coefficients  as  Eq. 
(23),  where  axis-1  denotes  the  coordinate  nonnal  to  the  plane  of  isotropy  [2]. 


cr, 


<Jn 


cr. 


cr,, 


<J< 


a, 


6j 


Qn  Qn  023  0  0  0 

Q22  023  0  0  0 

022  0  0  0 

1(011-012)  0  0 

066  0 

sym.  0r,6 


(23) 


As  the  lamina  we  are  dealing  with  can  be  treated  as  two  dimensional,  we  introduce  the 
plane  stress  state  assumption  (cr33  =0).  We  express  the  stress-strain  relationship  with  two 

separate  equations,  one  for  in-plane  components  and  the  other  for  transverse  ones.  In¬ 
plane  stress-strain  relations  in  on-axis  coordinates  (1 -,2-,3-axis)  are 


V 

V 

I 

LO 

LO 

to 

O 

_ 1 

£\ 

£\ ' 

°2 

>  =  < 

cr2 

.  = 

022  0 

•< 

£2 

=  0, 

£2 

CT6. 

,ri2. 

_sym.  066 

7l2. 

7X2 , 

(24a) 


where,  the  0,/s  are  expressed  in  terms  of  engineering  constants  in  the  following  manner 


0n  = 


1  ~  V12V21 


022  = 


1  ~  V12V21 


_  ^12^2  _  S'* 

&12  ,  9  zl-66  ^1: 

l-W2i 


(24b) 


The  transverse  shear  stress-strain  relationship  is 


a4 


23, 


066 

0 


0 


(25a) 
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where,  the  Qy  s  are  expressed  in  engineering  constants  as 


066 


0,2.  \(Qx~Qn) 


]_ 

2 


E\  ^2 
1  -^12^2, 


(25b) 


There  are  four  independent  material  constants  (EvE2,vn,Gn  )  for  the  case  of  thin  lamina 
assumed  in  the  plane  stress  state,  as  can  be  seen  in  Eqs.  (24)  and  (25). 

In  most  structural  application  of  composite  structures,  the  orthotropic  laminae  are 
stacked  into  a  laminate  with  a  certain  rotation  (or  lamination)  angle.  Fig.  3  shows  the 
relationship  between  the  material  axis  ( 1-2  axis,  on-axis)  and  the  laminate  axis  (x-y  axis, 
off-axis). 


Fig.  3  Positive  rotation  between  material  1-2  axis  (on-axis)  and  x-y  axis 


In-plane  stress  transformation  between  the  1-2  axis  and  x-y  axis  is  described  by 


cos2  9 

sin2  9 

2  cos  9  sin  9 

°x 

°x 

0-2 

>  = 

sin2  9 

cos2  9 

-2  cos  9  sin  9 

< 

>  =  T,< 

°y 

r!2. 

-cos  9  sin  9 

cos  9  sin  9 

cos2  9 -sin2  9 

.V 

(26) 


For  the  strain  transformation,  we  must  be  attentive  to  the  difference  of  strain  tensor  and 
engineering  strain.  The  strain  transformation  is  in  Eq.  (27). 


< 


cos2  9 

sin2  9 

2  cos  9  sin  9 

£x  ' 

sin2  9 

cos2  9 

-2  cos  9  sin  9 

< 

=  Tr 

£y 

-cos  <9  sin  <9 

cos  9  sin  9 

cos2  9 -sin2  9 

l  /2j 

r*y/ 

l  J 

(27) 


If  we  introduce  the  scaling  matrix,  R, 
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(28) 


R  = 


1  0  0 
0  1  0 
0  0  2 


then. 
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(29) 


The  off-axis  stress-strain  relation  for  in-plane  components  are  expressed  as  follows: 
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(30) 


where,  Q,  is  the  in-plane  off-axis  lamina  stiffness  matrix  and  its  elements  are  calculated 
from  the  rotation  angle  and  on-axis  stiffness,  Qfs. 

Qu=Qn  cos4  9  +  2(012  +  2066) sin2  9 cos2  9  +  022  sin4  9 
Qu  =  (0n  +Qn  -4066)  sin2  <9  cos2  9  +  012  (sin4  9  +  cos4  9) 

Q22=Qusin49  +  2(Ql2+2QJsin29cos29  +  Q22cos49 
Qu,  =  (Qn~  Qn  -  2066)  sin  6  cos3  9  +  (012  -  Q22  +  2066)  sin3  9  cos  9 
Q26  =  (0ii  -  0i2  -  2066)  sin3  9  cos  9  +  (012  -  022  +  2066)  sin  9  cos3  9 
066  =  (0ii  +  022  -  202  -  2066)  sin2  9  cos2  9  +  066  (sin4  9  +  cos4  9) 

Off-axis  transformations  of  transverse  shear  stress  and  strain  are  in  Eq.  (32)  and 
(33),  respectively. 
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(32) 
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(33) 


Similar  to  the  in-plane  stress-strain  relationship,  the  off-axis  transverse  stresses  and 
strains  are  related  as 


I  T* 

I  T„ 


=  T< 


-1  *"l3 


=  Ts-lQsV13\  =  Ts-lQsTs 


23 , 


rx: 

I  Yv. 


|  rxz 
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7l3 
1/23, 

055  054 

sym.  Q44 

where,  Qs  is  the  transverse  shear  off-axis  lamina  stiffness  matrix  and  its  elements  are 
calculated  from  the  rotation  angle  and  on-axis  stiffness,  Qs  s. 


055  =  066  cos2  #  + 1  (Q  j  -  Qv_ )  sin2  # 

054  =-066  COS# sin #  + 1(0!  1-0,2 ) cos# sin#  (35) 

044  =066  Sin2#  +  ^(0u-£>12)cos2# 


Once  several  laminae  are  stacked  to  construct  a  laminate  as  shown  in  Fig  (4),  we 
can  calculate  the  moment  resultant  by  summing  the  stresses  in  each  lamina.  Deformation 
represented  by  the  curvature  (  kx,  Ky,  kxv  )  is  related  to  the  moment  resultant 

( Mx , M  ,  Mxv  )  through  the  bending  stiffness  of  laminate,  DB ,  as  in  Eq.  (36). 
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In  a  similar  way,  the  transverse  shear  force  resultant  ( Qx,,Qyz )  is  expressed  in  terms  of 
the  transverse  shear  strain  ( yxz ,  y  )  in  the  laminate  with  the  transverse  shear  stiffness,  Ds  , 
which  is  shown  in  Eq.  (37). 


je„ 

[Q* 


=  TQs I  (z 


k=l  k 


k+l 


(37) 
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Fig.  4  z-coordinate  of  Lamina  in  a  Laminate 

D.  EIGENVALUE  ANALYSIS  AND  SENSITIVITY  OF  EIGENVALUE 

Assuming  that  the  displacement  response  is  harmonic,  the  equation  of  motion  in 
Eq.  (6)  can  be  expressed  as  Eq.  (38),  the  so  called  structural  eigenproblem. 

0}  (38) 

where,  A .  and  (f  are  y'-th  eigenvalue  and  eigenvector,  respectively.  The  eigenvector  is 

nonnalized  with  respect  to  mass  matrix.  The  mass  and  stiffness  matrices  in  Eq.  (38)  are 
obtained  by  assembling  the  element  matrices  in  Eqs.  (16)  and  (22),  and  boundary 
conditions  are  applied.  The  lowest  five  eigenvalues  and  corresponding  eigenvectors  are 


16 


calculated.  They  are  used  to  construct  the  perfonnance  function  for  optimization  together 
with  the  experimental  natural  frequencies  and  mode  shapes. 

As  the  four  mechanical  constants,  El,E2,vu,Gl2,  are  treated  as  design  variables 
for  optimization,  it  is  meaningful  to  investigate  the  sensitivity  of  natural  frequencies  with 
respect  to  these  design  variables  [15].  Differentiating  the  eigensystem,  Eq.  (38),  with 
respect  to  a  design  variable,  6 ,  yields 


dK 

~dd 


<t>,+K 


d</>j 

86 


8A,. 

— -Md),  +  A, 

86  1  1 


8M  8(f> 
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(39) 


Pre-multiplying  (j)J  to  above  Eq.  (39)  gives  us 
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8A , 
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86 
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(40) 


If  we  note  that  the  relationships  (  (j)J =  1  ;  K  and  M  are  symmetric;  and 
( K  -  AjM)^j  =  0  ),  Eq.  (40)  reduces  to 
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=  6.  — 6. 

’  86  1 


(41) 


for  the  case  of 


8M 

~86 


0, 


The  normalized  sensitivity  of  the  eigenvalue  is  defined  as  Eq.  (42): 

8Aj  _  l  f  8Aj ' 

86  Aj  k86  j 


(42) 


We  are  interested  in  the  normalized  eigenvalues  and  its  sensitivity  with  respect  to  the 
nonnalized  design  variable,  6  =6  /  60. 


8Aj 


8Aj  86 
86  86 


8  A.. 


86 


J00=- 


A, 


8K 


6.  - — 6. 

V  7  86  1 J 


6n 


(43) 


Eqs.  (16)~(18)  gives  us 
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8K‘ 
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The  normalized  eigenvalue  sensitivity  with  respect  to  the  nonnalized  design  variables 

depends  on  the  eigenvectors  and  the  derivatives  of  stiffness  matrices,  — and  L . 

89  89 

According  to  the  description  for  the  laminate  stiffness,  we  can  see  that  derivatives  of 
stiffness  matrices  are  expressed  eventually  in  tenn  of  derivatives  of  the  on-axis  stiffness, 
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III.  OPTIMIZATION 


A.  OVERALL  PROCEDURE 

Overall  procedure  for  obtaining  the  mechanical  property  constants  is  shown  in  the 
Fig.  5.  Four  mechanical  constants  are  treated  as  design  variables  for  optimization. 
Arbitrary  initial  values  are  assumed  for  the  design  variables.  Modal  parameters,  which 
are  the  natural  frequencies  and  mode  shapes  in  this  study,  are  calculated  with  these  initial 
values  for  the  system  model  of  the  specimen. 


Fig.  5  Schematic  procedure  to  obtain  mechanical  properties  through  optimization 
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These  modal  parameters  are  combined  with  the  experimental  ones  to  yield  the 
performance  index  to  be  minimized.  Until  satisfactory  minimization  of  performance 
index  is  obtained,  the  design  variables  are  updated  repeatedly.  Natural  frequencies  and 
mode  shapes  are  calculated  at  every  iteration.  Once  optimization  reaches  the  goal,  the 
natural  frequencies  of  the  mathematical  model  are  close  enough  to  the  experimental 
results.  The  design  variables  at  this  step  are  taken  as  the  desired  mechanical  properties  of 
the  lamina. 

It  is  required  to  have  experimental  modal  parameters  to  make  the  performance 
function  for  optimization.  At  present,  this  study  is  primarily  concerned  with  building  a 
computational  tool  for  mechanical  property  identification.  The  experimental  data  are 
simulated,  and  were  generated  from  the  analysis  results  with  reference  mechanical 
properties.  This  step  shall  be  replaced  when  experimental  data  are  available. 


B.  SPECIMEN  DESCRIPTION 

In  this  study,  the  system  model  is  a  laminated  plate  in  cantilever  configuration.  Its 
planform  is  as  shown  in  Fig.  6.  The  plate  is  modeled  with  a  nine-node  isoparametric 
element  based  on  the  first  order  shear  deformable  Reissner-Mindlin  plate  theory.  The 
laminate  will  be  constructed  with  orthotropic  lamina  whose  mechanical  properties  are  to 
be  identified.  The  lay-up  angle  of  each  lamina  are  specified  and  the  density  and  ply 
thickness  are  assumed  known.  Table  1  shows  the  reference  properties  for  graphite/epoxy 
lamina.  These  properties  would  be  obtained  through  a  serious  of  static  coupon  tests. 


Table  1.  Reference  mechanical  Proper  ties  of  Lamina 


E, 

122.5 

GPa 

E, 

7.929 

GPa 

U? 

0.329 

- 

3.585 

GPa 

Thickness  per  ply 

0.15 

mm 

Density  of  lamina 

1500 

kg/m 3 
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Clamped  Boundary 


Fig.  6  Description  of  the  Specimen 


C.  PERFORMANCE  FUNCTION  AND  OPTIMIZATION 


The  analytical  modal  parameters  are  used,  together  with  the  experimental  modal 
parameters,  to  construct  the  performance  function  to  be  minimized.  Design  variables  are 
updated  until  the  perfonnance  index  becomes  sufficiently  small.  Optimization  scheme, 
‘ fmincon  in  MATLAB’,  uses  a  sequential  quadratic  programming  method  and  BFGS 
formula  to  update  an  estimate  of  the  Hessian.  Details  of  the  scheme  can  be  found  in 
references  [16,  17]. 

Restrictions  on  the  engineering  constants  (Eq.  49)  which  come  from 
elastomechanical  constraints  are  handled  as  inequality  constraints  in  fmincon. 
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It  is  important  to  have  a  physically  reasonable  performance  function  for  this  kind 
of  problem.  The  mode  shape  for  ith  mode  obtained  experimentally  is  denoted  as  <f>eu)  and 

corresponding  analytical  one  is  (f>aU)  .  The  closeness  of  these  two  vectors  can  be 
represented  by  the  angle,  6j ,  between  them.  It  is  known  as  the  MAC ,. 


MAC,  = 


(0.)(ArA) 


=  cos2  Gj 


(50) 


In  the  case  that  a  specific  mode  shape  calculated  is  quite  different  from  one  obtained  in 
experiment,  it  is  of  no  use  to  try  to  match  the  natural  frequencies  from  analysis  and 
experiment.  Therefore  the  frequency  differences  for  each  mode  are  weighted  with  MAC s. 


J^MAC,)- 

i= 1 


I7.<',-/.(0Y 


/. 


(0 


(51) 


where,  faU)  and  fe(,)  are  natural  frequencies  from  analysis  and  experiment,  respectively. 
Nm  is  number  of  modes  in  consideration  and  five  in  this  study. 


D.  NUMERICAL  EXAMPLE 


As  a  numerical  example,  optimization  procedure  and  result  are  demonstrated  for  a 
laminate  whose  stacking  sequence  is  [±45°/ 0°2  /90°]s .  The  specimen  configuration  is  as 


in  Fig.  6.  Reference  properties  are  as  in  Table  1.  Modal  data  calculated  with  these 
properties  are  used  as  simulated  experimental  data.  A  starting  vector  of  mechanical 
constants  is  chosen  arbitrarily.  In  this  example,  it  is  given  as  follows: 


n-DV)=  4;  El 


12 


12 


A 


T7  R  ’ r  ’  r  R 

h2  vl2  tz[2 


=  [1.3;  1.3;  1.3;  1.3] 


where,  superscript  R  denotes  reference  value  as  in  Table  1.  Lower  and  upper  bounds  of 
the  nonnalized  design  variables  were  0.5  and  1.5,  respectively.  Table  2  shows  natural 
frequencies  with  optimized  mechanical  properties  together  with  the  reference  natural 
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frequencies.  Here  the  reference  natural  frequencies  are  treated  as  experimental  ones. 
Changes  in  design  variables  in  the  course  of  optimization  are  shown  in  Fig.  7  and  the 
performance  index  history  is  in  Fig.  8. 

Optimization  leads  the  natural  frequencies  of  the  specimen  to  the  target  value 
satisfactorily.  Table  2  shows  that  the  natural  frequency  for  each  mode  reaches  the 
corresponding  experimental  value.  As  we  do  not  have  in  this  study  any  error  or  noise  that 
usually  exists  in  real  experiments,  we  obtained  good  results  which  are  almost  the  same  as 
the  simulated  numerical  frequencies.  Fig.  9  shows  the  mode  shapes  obtained  with  the 
optimized  DV’s. 


Table  2.  Comparison  of  Natural  Frequencies  before  and  after  Optimization  (Hz) 


Mod 

With  Starting  DV’s 

With  Optimized  DV’s 

Experiment  (Simulated) 

1st 

70.12(114.4%) 

61.32  (100%) 

2nd 

270.57(113.8%) 

237.86  (100%) 

3rd 

423.30(114.2%) 

-► 

370.71  (100%) 

4th 

866.43  (114.0%) 

-> 

759.91  (100%) 

5th 

939.54(114.1%) 

— 

823.18  (100%) 
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Progress  of  Optimization 


Fig.  7  Optimization  progress,  Design  Variables,  starting  n_DV=  1.3*[1;  1;  1;  1] 
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Fig.  8  Optimization  progress,  Performance  Index,  starting  «_DF=1.3*[1;  1;  1;  1] 
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Progress  of  Optimization 
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Fig.  9a  Mode  shape  of  lsl  mode  with  reference  properties 
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Mode  2  (237.86  Hz) 


Mode  2  (237.86  Hz) 


Fig.  9b  Mode  shape  of  2nd  mode  with  reference  properties 
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Mode  3  (370.71  Hz) 


Mode  3  (370.71  Hz) 


Fig.  9c  Mode  shape  of  3rd  mode  with  reference  properties 
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Mode  4  (759.91  Hz) 


Mode  4  (759.91  Hz) 


0.15 


Fig.  9d  Mode  shape  of  4th  mode  with  reference  properties 
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Mode  5  (823.18  Hz) 


Mode  5  (823.18  Hz) 


Fig.  9e  Mode  shape  of  5th  mode  with  reference  properties 
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To  see  the  situation  with  different  starting  points  during  optimization,  two 
additional  starting  points  are  chosen.  Figs.  10a  and  10b  show  the  progress  of  optimization 
in  terms  of  design  variables.  For  the  case  of  Fig.  10a,  the  starting  vector  of  design 
variables  is  [0.7;  0.7;  0.7;  0.7],  i.e.,  70%  of  the  reference  values.  Optimization  yields  the 
values  of  [0.9999;  1.0005;  0.9982;  1.0006].  If  we  start  with  the  vector,  [  1 .4;  0.7;  0.6;  1.2], 
we  reach  the  result,  [1.0000;  0.9998;  1.0011;  0.9997].  It  cannot  be  mentioned  in  general 
that  any  starting  vector  may  give  the  desired  result,  [1.0000;  1.0000;  1.0000;  1.0000]. 
Concerning  the  initial  estimate  of  the  mechanical  properties,  which  is  necessary  for  the 
approach  in  this  study,  we  have  adequate  freedom  as  needed  to  choose  the  starting  vector. 

The  weighting  factor,  MAC,  in  the  performance  function  plays  a  role  to  reduce  the 
effect  of  specific  mode  of  which  the  mode  shape  from  the  mathematical  model  differs 
from  that  of  the  experiments.  The  experimental  mode  shape  vector  in  this  study  is 
constructed  with  the  z-displacement  at  the  center  node  of  finite  elements.  As  mention 
before,  it  is  obtained  not  from  experiment  but  numerically.  Table  3  shows  the  starting  and 
optimized  vector  of  design  variables  and  the  MAC s  for  the  mode  shape  with  several 
starting  vectors.  They  are  all  compared  with  the  mode  shapes  with  reference  properties 
which  are  treated  as  experimental  ones  in  this  study. 


Table  3  MAC s  for  mode  shapes  with  different  starting  vectors 


Case  1 

Case  2 

Case  3 

Starting  vector 

1 

Optimized  vector 

< 

1.3 

1.3 

1.3 

1.3 

>  -» < 

'o  o  o  o ' 

V. _ J 

| 

0.9999' 

1.0005 
0.9982  ’ 
1.0006 

1.4 

0.7 

0.6 

1.2 

1.0000' 

0.9998 

1.0011 

0.9997 

MAC 

of 

Starting  vector 

< 

1.0 

1.0 

1.0 

0.9996 

0.9995 

IMWI 

1.0 

1.0 

1.0 

0.9995 

0.9995 

> 

< 

1.0 

1.0 

1.0 

0.9827 

0.9818 

> 
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Normalized  design  variables 


Progress  of  Optimization 


Fig.  10a  Optimization  progress,  Design  Variables,  starting  w_DF=0.7*[l;  1;  1;  1] 
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Fig.  10b  Optimization  progress,  Design  Variables,  starting  n_DV=[\A;  0.7;  0.6;  1.2] 
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This  study  proposes  a  methodology  for  identifying  mechanical  properties  of 
orthotropic  lamina.  It  is  therefore  of  use  to  investigate  the  sensitivity  of  eigenvalues  or 
natural  frequencies  of  the  specimen  with  respect  to  the  mechanical  properties.  Table  4 
shows  them  using  Eq.  (43)  through  (48).  As  the  laminate  in  this  example  is  constructed 
with  many  layers  and  the  lamination  angles  are  ±45°,  0°,  and  90°,  the  stiffness  in  fiber 
direction,  E\,  has  prevailing  influence  on  all  modes  in  consideration.  Sensitivity  over  0.1 
is  shaded  in  the  table.  Gn  has  the  sensitivity  over  0. 1  only  for  1st  and  3ld  modes. 


Table  4.  Eigenvalue  sensitivities,  d^j/dO 


1st  mode 

2nd  mode 

3rd  mode 

4th  mode 

5th  mode 

Ei 

0.8589 

0.9221 

0.8338 

0.8627 

0.8436 

e2 

0.0403 

0.0393 

0.0327 

0.0716 

0.0633 

Vn 

0.0175 

-0.0174 

0.0083 

-0.0015 

0.0049 

Gn 

0.1008 

0.0386 

0.1334 

0.0656 

0.0931 

With  the  aid  of  Eq.  (52),  we  can  express  the  percent  of  changes  in  natural  frequencies 
according  to  the  changes  of  design  variables. 

L™=\\  +  *i=.(M)  (52) 

Lu  V  dd 

Table  5  shows  the  percent  changes  of  natural  frequencies  when  each  design  variable 
changes  10%.  If  we  change  Ej  by  10%,  the  natural  frequencies  up  to  5th  mode  changes 
by  approximately  4%.  Other  design  variables,  however,  seem  not  to  have  noticeable 
effect  on  the  natural  frequencies  compare  to  Ej. 


Table  5.  Changes  of  natural  frequencies  followed  by  10%  change  of  DV’s 


1st  mode 

2nd  mode 

3rd  mode 

4  th  mode 

5th  mode 

AE1=10% 

4.21  % 

4.51  % 

4.09  % 

4.22  % 

4.13  % 

AE2=10% 

0.20  % 

0.20  % 

0.16% 

0.36  % 

0.32  % 

A  Up  =10% 

0.09  % 

-0.09  % 

0.04  % 

-0.01  % 

0.02  % 

AG12=10% 

0.50  % 

0.19% 

0.66  % 

0.33  % 

0.46  % 
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The  major  influence  of  Ej  on  the  natural  frequencies  is  mentioned  for  the  laminate 
of  the  previous  example.  For  a  laminate  which  is  constructed  with  laminae  in  special 
orientations,  the  situation  becomes  different.  Table  6  and  7  show  the  eigenvalue 
sensitivity  for  the  laminates,  [±45°]2V  and  [0°]gr .  As  the  lamination  angles  are  tailored  for 

specific  stiffness,  we  can  see  major  influence  of  design  variables  other  than  E]  on  the 
eigenvalues  in  the  Tables  below.  Again,  sensitivities  whose  magnitudes  are  over  0.1  are 
shaded  in  the  tables. 


Table  6.  Eigenvalue  sensitivities,  dXj  jdd  for  laminate,  [±45°]2v 


1 st  mode 

2nd  mode 

3rd  mode 

4th  mode 

5th  mode 

Ei 

0.5050 

0.6579 

0.4767 

0.5233 

0.6485 

e2 

0.0534 

0.0532 

0.0415 

0.0482 

0.0667 

*12 

0.0144 

-0.0089 

0.0063 

0.0065 

-0.0036 

Gl2 

0.2130 

0.0348 

0.2530 

0.1947 

0.0378 

Table  7.  Eigenvalue  sensitivities,  dX j/dO  for  laminate,  [0°]gr 


1st  mode 

2nd  mode 

3rd  mode 

4th  mode 

5th  mode 

Ei 

0.7128 

0.3891 

0.0870 

0.7076 

0.6047 

e2 

0.0013 

0.0158 

0.8721 

0.0015 

0.0184 

*12 

0.0033 

0.0026 

0.0031 

0.0034 

0.0022 

G 12 

0.0010 

0.3702 

0.2232 

0.0069 

0.1172 
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IV.  DESCRIPTION  OF  COMPUTATIONAL  ROUTINES 


All  the  procedures  in  this  study  are  written  in  MATLAB.  Brief  descriptions  of  the 
major  steps  and  variables  are  given  in  this  chapter  and  source  listing  is  attached  in 
Appendix. 


A.  Main  Routine 

The  main  routine  (main _procedure  Jor  oplim.m)  defines  the  starting  vector  of  the 
nonnalized  design  variables  ( n_DV)  and  information  of  the  specimen:  specimen  length  in 
x-direction  (La)  and  in  y-direction  (Lb),  number  of  finite  element  elements  in  x-direction 
(Na)  and  in  y-direction  (Nb),  number  of  layer  in  the  laminate  (NT),  lamination  angles 
(Ang),  density  of  material  (rho),  and  unit  thickness  of  lamina  (tlayer).  Using  the 
information  of  the  specimen,  the  main  routine  calls  ‘ modeling  comp/at e. m ’  to  generate 
information  for  finite  elements  such  as  the  coordinate  of  the  grid  points  and  the  element 
connection  of  the  finite  elements,  etc. 

Experimental  results  such  as  natural  frequencies  and  mode  shapes  are  to  be 
imported.  As  explained  before,  reference  analytical  data  are  prepared  to  substitute 
temporarily  for  the  experimental  natural  frequencies  and  mode  shapes.  Upper  and  lower 
bounds  for  design  variables  (ub  and  lb)  and  several  option  parameters  (TolFun,  To/Con, 
and  etc.)  are  specified.  Constrained  minimization  routine,  fmincom.m,  is  called.  As  input 
parameters  for  this  routine,  two  additional  routines  are  prepared.  One  is  for  the 
performance  index  calculation  (pi Jv.m)  and  the  other  is  a  routine  to  define  nonlinear 
constraints  between  design  variables  (restr_eng_const.m). 


B.  Perforamnce  index  calculation 

The  routine,  pi Jv.m,  is  repeatedly  called  during  the  minimization  process  of  the 
performance  index.  This  routine,  pi Jv.m,  calls  the  routine,  sol Jv.m,  which  returns  the 
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analytical  natural  frequencies  and  the  mode  shapes.  The  performance  index  is  calculated 
using  the  experimental  and  analytic  modal  data.  The  variable,  hist,  keeps  the  record  of 
progress  of  optimization  such  as  performance  index  and  design  variables. 


C.  Solve  Eigenproblem 

The  routine,  sol Jv.m,  performs  element  matrix  generation,  and  assembly  to 
generate  global  matrices,  apply  boundary  condition,  and  call  eigenproblem  solver.  To 
find  lower  several  eigenvalues  and  mass  normalized  eigenvectors,  the  routine,  eigsn.m,  is 
are  used.  It  is  a  modified  routine  of  eig.ra  and  eigs.m.  Analytical  eigenvalues  and 
eigenvectors  up  to  a  certain  mode  (here  5th)  are  returned  to  the  routine,/;/ _Jv.m.  From  the 
eigenvector,  vertical  component  at  the  center  node  of  finite  elements  are  chosen  and 
rearranged  to  be  used  for  comparison  with  the  experimental  mode  shapes. 


D.  Finite  element  generation 

There  are  several  routines  for  finite  element  generation  as  follows: 

.  MeKe.m  generate  element  stiffness  and  mass  matrices 

.  dKedDVi.m  generate  derivatives  of  element  stiffness  matrix 

.  D  matrix.m  calculate  bending  stiffness  of  laminated  plates 

.  D_sen_matrix_wrt_DVi.m  calculate  derivatives  of  D  matrix 

.  sen_eigvalue.m  calculate  sensitivities  of  eigenvalue 

.  shape _iso9.m  calculate  value  and  derivatives  of  the  shape  function 

of  9-node  isoparametric  element 
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V.  CONCLUSION  AND  RECOMMENDATION 


A  method  of  obtaining  the  mechanical  properties  of  the  orthotropic  lamina  is 
presented  along  with  computational  routines.  Differences  between  the  natural  frequencies 
from  mathematical  model  and  those  from  experiment  are  minimized  by  updating  the  four 
mechanical  stiffness  of  lamina.  The  frequency  difference  in  each  mode  is  weighted  based 
on  the  modal  assurance  criteria.  A  simple  vibration  test  to  obtain  the  natural  frequencies 
and  mode  shapes  can  be  a  substitute  for  a  series  of  static  coupon  test  to  characterize  the 
four  mechanical  stiffness  constants. 

This  study  utilizes  finite  element  analysis  using  nine-node  Reisnner-Mindlin  plate 
element,  the  classical  lamination  theory,  and  the  first  order  shear  deformation  theory. 
Each  procedure  is  coded  in  MATLAB,  which  is  included  in  this  report.  The  MATLAB 
built-in  function,  fmincon,  is  used  to  minimize  the  performance  index.  A  numerical 
example  is  given  to  demonstrate  the  performance  and  usefulness  of  this  scheme. 

It  is  an  inverse  problem  to  find  the  four  design  variables  which  can  match  the 
dynamic  response  of  a  mathematical  model  of  a  specimen  with  the  real  experimental  one. 
It  is  necessary,  therefore,  to  have  experimental  natural  frequencies  and  mode  shapes  in 
addition  to  the  computational  tools.  Future  work  to  get  the  experimental  data  are 
recommended  to  complete  this  work  for  characterization  of  the  mechanical  stiffness 
constants  of  orthotropic  lamina. 
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APPENDIX 


A.  mainprocedureforoptim 


% -  main_procedure_for-optim.m  -- 

"o 

%  main  program  for  mechanical  property  identification 

"o 

%  for  laminate,  [+/-45/0/0/ 90 ] s 


%  Overall  procedure 


0)  define  the  design  variables:  DV's 

-  with  proper  nomalization 

1)  describe  the  specimen:  run  "modeling_complate . m" 

-  composite  laminate,  grid  point  definition,  element  connetion, , , 

-  declare  global  variable,  ***************** 

2)  import  experimental  results 

-  natural  frequencies  and  mode  shapes 

-  declare  global  variables:  f_exp,  v_exp 

3)  call  optimization  routine:  call  "fmincon.m"  and  "@pi_fv.m" 

4)  in  the  performance  index  function,  "pi_fv.m" 

4.0)  decleare  global  variables 

4.1)  element  matrices:  call  "MeKe.m" 

-  mass  matrix,  stiffness  matrix 

4.2)  assemble  of  element  matrices  and  apply  boundary  condition 

4.3)  solve  eigenproblem :  call  "eigsn.m" 

4.4)  extract  eigenvalue  and  eigenvector (selected  dof) 

4.5)  construct  the  performance  index 

-  compare  the  diffenece  in  experiment  and  analysis 

-  minimization  of  difference  in  eigenvalue  and  . . . 

angle  between  eigenvector  and  mode  shape 

8)  postprocess  the  resultes 

-  plot,  etc. 


%  Input  information 


[Laminate  comfiguration] 

-  La:  length  in  x-direction  (meter) 

-  Lb:  length  in  y-direction  (meter) 

-  Na:  #  of  element  in  x-direction 

-  Nb:  #  of  element  in  y-direction 

-  Nl:  #  of  layer  in  composite  laminate 

-  Ang:  lamination  angle 

-  rho :  material  density  ( kg/meter A3 ) 

-  tlayer:  lamina  thickness (meter) 

[Finite  element  model] 

-  [Coord-grid]:  grid  point  coordinates,  [xl,yl;  x2,y2;;;;  xn,yn] 

-  [Elem] :  element  connectivity, 

[EID_1 , 9-gid's; EID_2 ,9-gid's;;;;;  EID_n, 9-gid ' s ] 
e.g.(for  EID=A  shown  below) ~  Elem= [A, 1 , 15, 17, 3, 8, 16, 10, 2,  9] 

-  [Lam] :  laminate  information,  [EID,  layer  id,  thickness,  angle] 

e.g.  [Eid_l,  1_1,  t_l,  a_l;  Eid_l,  layer_2,  t_2,  a_2;;;; 

Eid_2,  1_1,  t_l,  a_l;  Eid_l,  layer_2,  t_2,  a_2; ;;;;;;] 

-  [Be] :  applied  boundary  condition,  fixed  @x=0 

[gid,  zeros  to  constrain  as  many  as  dof  per  grid) ~  BC 
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%  [Design  variables  =  initial  values  of  material  properties] 

%  -  [DV] :  material  properties:  El , E2, vl2, G12  (SI  unit) 

"o 

%  [Eexperimental  data] 

%  -  [f_exp] :  natural  frequencies  in  Hz  (1  x  nmode) 

%  -  [v_exp] :  mode  shape  (exp_dof  x  nmode) ,  exp_dof=Na*Nb  @center  node 

%  *  each  mode  shape  is  in  the  order  as  follows: 

%  [^center  of  EIDl,@center  of  EID2 , , , , , @center  of  EID(Na*Nb)]' 

-6 

%  [Analytical  data] 

%  -  [f  ana] :  eigenvalues  in  Hz  during  updating  design  variables 

%  -  [v_ana] :  eigenvectors  with  selected  dof ' s  (exp_dof  x  nmode) 

"6 

%  [Performance  index] 

%  -  J=sum (MAC (i ) *sqrt ( ( f_ana ( i) -f_exp (i ) ) A2/f_exp (i ) A2) ) 

%  -  where,  MAC ( i) = ( v_ana ( : , i ) ' *v_exp ( : , i ) ) A2/ . . . 

%  ( (v_ana ( : , i) ' *v_ana ( : , i) ) * (v_exp ( : , i) ' *v_exp ( : , i) ) ) 

"6 

%  17  Spetember  2007  (jkr) 

"6 

clear  all 
%  clear 

%  clear  global 
"6 

global  hist;  %  history  of  optimization,  fmincon.  accumulated  @  pi_fv.m 


%  define  design  variables  and  use  as  initial  values 


%  define  initial  value  of  design  variables,  DV 
"6 

global  DV_r; 

"6 

DV_r= [122 . 5e9;  7.929e9;  0.329;  3.585e9];  %  [El  E2  vl2  G12] 

"6 

n_DV=0.7* [ 1 ;  1 ;  1 ;  1  ]  ; 

%  n_DV= [ 1 . 4 ; 0 . 7 ; 0 . 6 ; 1 . 2 ]  ;  %  arbitrary  initial  value 
%  n_DV=l . 3* [ 1 ; 1 ; 1 ; 1 ] ; 

"o 

%  material  property  DB 


"o 

mid 

rho 

El 

E2 

vl2 

G12 

"o 

0 

1500 

122. 5e9 

7 . 92  9e9 

0.329 

3 . 585e9 ; 

%  t=0. 150mm  Baseline 

-6 

1 

1600 

207e9 

5e9 

0.25 

2 . 6e9 

%  Gr/Ep  (Jones) 

"o 

2 

1600 

122 . 4552e9 

7 . 92  92  5e9 

0.329 

3.5854e9; 

%  t=0. 125mm  Ref. 

"o 

3 

1450 

122 . 4552e9 

7 . 92  92  5e9 

0.329 

3.5854e9; 

%  t=0.13mm,  jkr 

"o 

4 

1500 

108. 4e9 

7 . 703e9  0. 

.3193  2. 

.  776e9;  % 

t=0 . 1 5mm,  danbi 

"o 


%  description  of  specimen,  dimensions,  lay-up,  and  etc. 

"O 

-6 


La=0 .15; 

Lb=0 . 1; 

Na=6  ; 

Nb=4  ; 

Nl=10 ; 

Ang= [ 45  -45  0  0  90  90  0  0  -45  45] ; 

rho=1500; 

tlayer=0.150e-3; 


%  length  in  x-direction  (meter) 
%  length  in  y-direction  (meter) 
%  #  of  element  in  x-direction 
%  #  of  element  in  y-direction 
%  #  of  layers  including 
%  lamination  angle 
%  density  (kg/mA3) 

%  lamina  thickness  (meter) 
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%  finite  element  model  data  using  9-node  isoparametric  element 


global  Coord_grid  Elem  Lam  Be; 


%  modeling  results 


[Coord_gr id, Elem, Lam, Be ] =modeling_complate (La, Lb, Na, Nb, N1 , Ang, t layer) ; 

-6 

%  parameters  for  finite  element  analysis  model 


nnel=9; 

ndof=3; 

edof=nnel*ndof  ; 
ngr id= (2*Na+l)* (2*Nb+l) ; 
t_dof =ndof *ngrid; 
n  elem=Na*Nb; 


%  number  of  nodes  per  element 
%  dof's  per  node,  [phil , phi2 , w] ' 

%  dof's  per  element 

%  total  number  of  nodes  (grid  point) 
%  total  dof 
%  #  of  elements 


global  s_info  m_info; 

-6 

s_inf o= [La, Lb, Na, Nb, Nl, Ang, rho, tlayer ] ;  %  specimen  information,  s_info 

m  info= [nnel, ndof , edof , ngrid, t  dof,n  elem];  %  model  information,  m_info 


%  import  experiment  results 


%  As  there  are  no  available  experimental  modal  data,  fictitious 
%  experimental  data  are  calculated  with  the  reference  properties 
-6 

%  -  for  [+45/-45/0/0/90 ] s  laminate 

%  in  the  file:  "ref_laml_AR150 .mat" 

%  data:  f_r  v_r  mshp_gid_r  lambda_r  Phi_r  dn_lamdn_DV_r 

"o 

%  where,  f_r:  reference  natural  frequency  up  to  5th  mode 
%  v_r :  reference  mode  shape  up  to  5th  mode 

%  mshp_gid_r:  mode  shape  for  plotting 

%  lambda_r:  reference  eigenvlaues 

%  Phi  r:  reference  eigenvectors 

%  dn_lamdn_DV_r :  eigenvalue  sensitivity  @ref  properties 

-6 

global  f  exp  v  exp  mac;  %  experimental  results,  weighting 
global  nmode; 


nmode=5 ; 


%  number  of  modes  being  considered 


load  ref_laml_ARl 50  f_r  v_r  mshp_gid_r  lambda_r  Phi_r  dn_lamdn_DV_r ; 

"o 

f  exp=f  r;  v  exp=v  r; 


global  f_ana  v_ana  mshp_gid  lambda  Phi;  %  defined  in  the  routine,  pi_fv_r2 


%  for  optimization  using  fmincon 

%  -  performance  function,  pi_fv.m 

%  -  constraints  function,  restr_eng_const .m 

"6 

%  -  n_o :  optimized  design  variables,  normalized 

%  -  J_o :  minimized  performance  index 


lb=0 . 5* [ 1 ; 1 ; 1 ; 1 ] ; 
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ub-1.5* [ 1 ; 1 ; 1 ; 1 ] ; 

"6 

%  my_opt=optimset . . . 

%  ('Display', 'notify',' FunVal Check ' , 'on', 'TolFun',le-8, ’TolCon',le-8) ; 

my_opt=optimset . . . 

('Display', 'iter', ' FunVal Check ', 'on', 'TolFun',le-8, 'TolCon',le-8) ; 

"6 


dfile=input (' diary  file  name:  ',  's'); 

save(dfile,  'hist')  %  save  raw  data  in  'dfile.mat' 
diary  (dfile) 

disp  (  ' - ' ) 

disp (' optimization  using  [ 4 5/-4 5/0/0/ 90 ] s  laminate') 

disp  (  ' - ' ) 

disp( ['start  design  variables  =  [  ',  num2 str (n_DV ' ) ,  '  ]']) 

disp  (  ' - ' ) 

disp( ['start  natural  freqency  =  [  ',  num2 str ( f_exp ' ) ,  '  ]']) 

disp  (  ' - ' ) 

[n_o, J_o, exitflag_o, output_o] =. . . 

fmincon (@pi_fv, n_DV, [ ] , [ ] , [ ] , [ ] , lb, ub, @restr_eng_const, my_opt) ; 

disp  (  ' - '  ) 

disp(['optim  design  variables  =  [  ',  num2str (n_o ' ) ,  '  ]']) 

disp  (  ' - ' ) 

disp(['optim  natural  freqency  =  [  ',  num2 str ( f_ana ' ) ,  '  ]']) 

disp  (  ' - ' ) 

f_ratio= ( f_ana . / f_exp) *100; 

disp(['optim  natural  freqency  =  [  ',  num2 str ( f_ratio ' ) ,  '  ]']) 

disp  (  ' - ' ) 

disp ( '  END  ' ) 

disp  (  ' - ' ) 


diary  off 
-6 
"o 

%  optimization  history  plot 

"o 

[fid, message] =fopen (dfile, ' r ') ;  %  open  data  file 
-6 

%  check  file  open 

"o 

if  ~isempty(message); 

disp (' error !!  !  To  find/open  diary  file') 
break 

end 

"o 

f rewind (fid) ; 

"o 

%  move  to  location  of  iteration  data 

"o 

while  1 

n_line=fgetl (fid)  ; 

[aaa,  line_length] =size (n_line) ; 
if ( line_length  <  25); 
n_line=fgetl (fid) ; 

end 

if  n_line(2:25)  ==  'Iter  F-count  f (x) ' ; 

break 

end 

end 

"o 

%  Initialize  storage: 
xx_iter= [ ] ; 
xx_fcount= [ ] ; 
xx  f = [ ] ; 
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%  read  iteration  starting  status 

"o 

n_line=fgetl (fid)  ; 

xx_iter=eval (n_line (1:5)); 

xx_f count =eval (n_line ( 6 : 12 )  )  ; 

xx_f=eval (n_line (13:25)  )  ; 

xx(l,:)  =  [xx  iter,  xx  fcount,  xx  f ]  ; 

"o 

%  read  iteration  history 

"o 

n_iter=l ; 

n_line=fgetl (fid) ; 
while  (eval (n_line (1:5) ) ==n_iter ) 
xx_iter=eval (n_line (1:5)) ; 
xx_f count =eval (n_line (6:12)  )  ; 
xx_f=eval (n^line (13:25)  )  ; 

xx(l+n  iter,:)=[xx  iter,  xx  fcount,  xx  f] ; 
n_iter=n_iter+l ; 
n_line=fgetl (fid) ; 

[aaa, ok] =str2num (n_line (1:5)); 
if  ok==0;  break;  end; 

end 

"o 

fclose ( fid) ; 

"o 

%  optimization  history,  opt_hist 

%  opt_hist= [ iter# ,  Function_call_count,  f_value,  DV_1~DV_4,  fl~f5] 

"o 

opt_hist= [xx ( : , 1) , xx ( : , 2) , hist (xx ( : , 2 ) , : ) ] ; 

"6 

figure 

plot (opt^hist ( : ,  (4:7))) 
xlabel (' Iteration  Number') 
ylabel ( 'Normalized  design  variables') 
title (' Progress  of  Optimization') 
legend ( 'El ' , 'E2 ' , 'vl2 ' , 'G12 ' ) 
grid  on 
-6 

figure 

plot (opt^hist ( : , 3) ) 
xlabel (' Iteration  Number') 
ylabel (' Per formance  Index') 
title (' Progress  of  Optimization') 
legend ( ' PI ' ) 
grid  on 

"5 

figure 

plot (opt^hist  ( : ,  ( 8  : 12 ) ) ) 
xlabel (' Iteration  Number') 
ylabel (' Natural  Frequency (Hz) ') 
title (' Progress  of  Optimization') 
legend ( ' f 1 ' , ' f 2  ' ,  ' f 3 ' , ' f 4 ' ,  '  f 5  ' ) 
grid  on 

"o 

% -  end  of  main_procedure_for-optim.m  - 

-6 
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B.  modelingcomplate 


modeling_complate  .m 


function  [Coord_grid, Elem, Lam, Be] . . . 

=modeling_complate (La, Lb, Na, Nb, N1 , Ang, t layer ) 

"o 

%  Input  preparation  for  composite  plate:  La  x  Lb  clamped  at  root  (at  x=0) 


[Input  parameters] 

-  La:  length  in  x-direction  (meter) 

-  Lb:  length  in  y-direction  (meter) 

-  Na:  #  of  element  in  x-direction 

-  Nb:  #  of  element  in  y-direction 

-  Nl:  #  of  layer  in  composite  laminate 

-  Ang:  lamination  angle  (note:  reverse) 

-  rho :  material  density  ( kg/meter A3 ) 

-  tlayer:  lamina  thickness (meter) 

-  DV:  material  properties:  El , E2, vl2, G12  (SI) 

[Output  parameters] 

-  Coord_grid: (x,  y) ~  coordinates  of  grid  points  in  sequential  order 

-  Elem: (EID,  grid_id's)~  element  connectivity 

e.g.(for  EID=A  shown  below) ~  Elem= [A, 1 , 15, 17, 3, 8, 16, 10, 2, 9] 

-  Lam  :  (EID,  layer  id,  thickness,  angle) ~  laminate  information 

-  Be  : (grid  id,  zeros  to  constrain  as  many  as  dof  per  grid) ~  BC 


--14 - 21 

I  I 

===C===2  0 
I  I 

—  12 - 19 

I  I 

===B===1 
I  I 

—  10 - 17 

I  I 

===A===1 6 
I  I 

—  _8 - 15 


I 

=1  = 

I 

-+- 

! 

==== I ==== 

I 

-+- 

I 

I 


1,2,3,,,  :  GID 

-I - + 

A, B, C, , ,  :  EID 

-I - + 

- >x 


%  17  September  2007  (jkr) 

"o 

Coord_grid=zeros ( (2*Na+l) * (2*Nb+l) ,2) ; 
for  i=l : 2  *Na+l 

for  j  =1 : 2 *Nb+l 

Coord_grid ( (i-1) * (2*Nb+l) +j , : )  =  [La/2/Na* (i-1)  ,  Lb/2/Nb* (j-1)  ]  ; 

end 

end 


Elem=zeros (Na*Nb, 10) ; 
for  i=l : Na 

for  j  =1 : Nb 

fst= ( j-1) *2+ (i-1) *2* (2*Nb+l) +1;  %  first  grid  id  for  each  element 
Eid= ( i-1 ) *Nb+ j ;  %  element  id 
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Elem (Eid, :)=... 


[Eid, ...  % 

fst,  fst+2* (2*Nb+l)  ,  fst+2* (2*Nb+2) ,  fst+2,...% 

fst+ (2*Nb+l) ,  fst+2* (2*Nb+l) +1,  fst+ (2*Nb+l) +2,  fst+l,...% 

fst+ (2*Nb+l) +1] ;  % 

end 

end 

clear  fst 

"o 

clear  Lam 
for  i=l : Na 

for  j  =1 : Nb 

L_elm= [ (i-1 ) *Nb+j ,  Nl,  tlayer,  Ang(Nl)];  %  top  layer  of 
for  k=2:Nl;  %  loop  to  stack  layer 

temp= [ (i-1) *Nb+j , (Nl+1) -k,  tlayer,  Ang ( (Nl+1 ) -k) ] ;  % 

%  Lam ( ( i-1 ) *Nb+ j : ( i-1 ) *Nb+ j , : ) = [Lam;  temp]; 

L_elm= [L_elm; temp] ;  %  stack  layer 
%  [Lam;  (i-l)*Nb+j,  (Nl+l)-k,  tlayer,  Ang ((Nl+1) - 

end 

Lam ( (Nl*Nb* (i-1) +N1* ( j-1) +1)  :  ( N 1 *  Nb  * (i-1) +j *N1) , : ) =L_elm; 

end 

end 

-6 

%  Boundary  Codition:  Clamped  along  the  edge,  x=0 

"o 

Bc=zeros (2*Nb+l , 4) ; 
for  i=l : 2  *Nb+l 

Be (i, : )=[i  0  0  0  ] ; 

end 

"o 

clear  Eid  L_elm  temp  i  j  k 

"o 

% - end  of  modeling_complate . m  - 


element  id 
corner  gid 
side  gid 
center  gid 


an  element 
next  layer 


k)  ]  ; 
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C.  pi_fv 


"o 

"o 

function  J_all=pi_fv (n_DV) 

"o 


%  performance  index  calculation 

"o 


pi_fv .m 


%  17  September  2007  (jkr) 


global 

global 

global 

global 

global 

global 

"5 

"5 


Coord_grid  Elem  Lam  Be;  %  model  data 

nmode ; 

f_exp  v_exp  mac;  %  experimental  results,  weighting  (from  main) 
f_ana  v_ana  mshp_gid  lambda  Phi;  %  to  be  used  in  main  routine 
DV_r;  %  reference  value  of  design  variable 

hist; 


%  obtain  eigenvalue,  eigenvector,  and  mode  shape 


[ f_ana, vana, mshp_gid, lambda, Phi ] =sol_fv (n_DV) ; 

"6 

"o 

%  performance  index  construction,  J=sum{MACi*del ( f i ) } 

"6 

-5 

for  n=l: nmode 

mac (n, 1 ) = (v_ana ( ; , n) ' *v_exp (;,n))A2/... 

( (v_ana ( ; , n) ' *v^ana ( : , n) ) * (v_exp ( : , n) ' *v_exp ( : , n) ) ) ; 
delf (n,  1 ) =sqrt ( (f_ana (n) -f^exp (n) ) A2/ f_exp (n) A2 ) ; 

end 

"o 

%  J_all=ones ( 1 , nmode) *delf;  %  without  mode  shape  weighting 
J_all=mac ' *delf ;  %  with  mode  shape  weighting 

-6 
"o 


%  convergence  history 

%  each  row  for  iteration  containg,  [PI, n_DV, f ' s ] 

"6 

"o 

c_hist= [ J_all, n_DV' , f_ana ' ] ; 
hist= [hist; c_hist] ; 

"o 

% - end  of  pi_fv.m - 


"o 

"o 


weighting 
del  freq 
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D.  restengconst 


% -  rest_eng_const .m  - 

"o 

function  [c, ceq] =restr_eng_const (x) 

"o 

global  DV_r  %  baseline  engineering  constants,  [El  E2  vl2 

"o 

%  engineering  constants: 

%  El=x ( 1 ) *DV_r ( 1 ) ,  E2=x (2) *DV_r (2) ,  vl2=x (3) *DV_r (3) , 

"O 

%  input  x:  normalized  engineering  constants 
%  vl2-sqrt (E1/E2)  <  0 

%  or  (vl2) A2-  (E1/E2 )  <  0 

%  or  x (2) /x (1) *x (3) A2-DV_r (1) /DV_r (2) /DV_r (3) A 

"5 

%  output  c:  inequality  constraints  for  fmincon,  c(x)<=  0 
%  ceq:  equality  constraints  for  fmincon,  ceq  =  0 

"o 

%  17  Spetember  2007  (jkr) 

-6 

c=x (2)/x(l)*x(3) A2-DV_r (1) /DV_r (2) /DV_r (3) A2; 
ceq= [ ]  ; 

"o 

% - end  of  rest_eng_const .m  - 


G12  ] 

G12=x (4) *DV_r (4) 

2  <  0 
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E, 


sol  fv 


% -  sol_fv.m  - 

"o 

function  [ f_ana, vana, mshp_gid, lambda, Phi ] =sol_fv (n_DV) 

"o 

%  obtain  eigenvalue  and  eigenvector  for  given  n_DV 

"o 

%  17  September  2007  (jkr) 

-6 

global  DV_r;  %  reference  mechnical  properties 
global  Coord_grid  Elem  Lam  Be;  %  model  data 

global  s_info  m_info; 
global  nmode; 

global  f  exp  v  exp  mac;  %  experimental  results,  weighting 


La=s_info ( 1 ) ;  Lb=s_info (2 ) ;  %  plate  size 

Na=s_info (3 ) ;  Nb=s_info ( 4 ) ;  %  #  of  elements  in  each  direction,  x  &  y 
Nl=s_info ( 5 ) ;  Ang=s__inf o ( 6 : 5+N1 ) ;  %  #  of  layer  and  lamination  angle 

rho=s_inf o ( 6+N1 ) ;  tlayer=s_info ( 7+N1 ) ;  %  density  and  lamina  thickness 


nnel=m_info ( 1 ) ; 
ndof=m_info (2)  ; 
edof=m_info (3)  ; 
ngrid=m_inf o ( 4 )  ; 
t_dof=m_info (5)  ; 
n  elem=m  info (6); 


%  number  of  nodes  per  element 
%  degrees  of  freedom  per  node 
%  degrees  of  freedom  per  element 
%  total  number  of  grids  (nodes) 

%  total  dof 
%  #  of  elements 


DV=n_DV . * DV_r ;  %  convert  normalized  values  to  physical  ones 

"6 

-6 

%  finite  element  matrices  constuction 

"O 

%  assemble  element  matrix 
"6 

M=zeros (t_dof , t_dof ) ;  %  initialization  of  M-matrix 

K=zeros (t_dof , t_dof ) ;  %  initialization  of  K-matrix 

-6 

for  Eid=l : n_elem; 

"6 

[Me, Ke] =MeKe (Eid, DV) ; 

"6 


gids=Elem (Eid, 2:10) ; 
idx=feeldof (gids, nnel, ndof ) 

-6 

M=feasmbll (M, Me, idx) ; 
K=feasmbll (K, Ke, idx) ; 


%  grid  id's  of  the  element 
%  system  dof's  of  the  element 

%  assemble  of  system  mass  matrix 
%  assemble  of  system  stiffness  matrix 


end 

"6 

%  apply  boundary  condition  (cantilever  plate  -  LHS  clamped) 

-6 

[nr ,  nc] =size (Be) ; 

nbc=nr* (nc-1 )  ;  %  total  #  of  dof  constrained 

Ma=M (nbc+1 : t_dof , nbc+1 : t_dof ) ;  %  patitioned  mass  matrix (BC's  applied) 

Ka=K (nbc+1 : t_dof , nbc+1 : t_dof ) ;  %  patitioned  stiffness  matrix (BC's  applied) 

"6 

clear  Eid  gids  nr  nc 

"6 

"6 
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%  eigenvalue  analysis 


[Phi_, D_] =eigsn (sparse (Ka) , sparse (Ma) , nmode) ; %  few  modes  with  normalization 
lambda_=diag (D_) ; 

for  ii=l: nmode;  %  re-order  increasing  eigenvalue 

lambda ( ii , 1 ) =lambda_ ( nmode- ii+1 ) ; 

Phi ( ; , ii) =Phi_ ( : , nmode-ii+1 ) ; 
end 


%  %  Alternative  routine  for  solving  eigen  problem  above 
%  [ lambda_, Phi_, Psi ] =eign (Ka, Ma) ; 

%  lambda=lambda_ ( 1 : nmode) ; 

%  Phi=Phi  (:, 1 : nmode); 


%  rearrange  eigenvector  to  plot  the  mode  shapes 

"O 

"o 

mshp_gid=zeros (ngrid, nmode) ;  %  initialize 

v_ana=zeros (n_elem, nmode) ;  %  initialize 

-6 

eigmtrx_=Phi (:, 1 : nmode) ;  %  eigenvectors  up  to  'nmode'  modes 

eigmtrx= [ zeros (nbc, nmode) ; eigmtrx_] ;  %  add  zeros  at  fixed  boundary  dof's 

"o 

%  construct  mode  shape  with  w's  (  get  rid  of  rotation  dof's  ) 

-6 

%  mode  shape  along  first  along  y  then  increase  x  (  along  gid  ) ,  mshp_gid 

"o 

for  i=l : ngrid 

mshp_gid ( i, : ) =eigmtrx (ndof *i , : ) ;  %%%  for  mode  shape  plot 

end 

"6 

"o 

%  arrange  analysis  results  to  compare  with  experiments 

"o 

"6 

f_ana=sqrt ( lambda ( 1 : nmode) ) /2 /pi ;  %%%  natural  frequencies  from  analysis 

"6 

%  eigen  vectors  @  center  of  element:  nelem= (Na) x (Nb ) 

"6 

for  i=l : Na 
for  j  =1 : Nb 

v_ana ( j +Nb* ( i-1 ) , : ) =mshp_gid ( 2* j + (2 *Nb+l ) * ( 2* i-1 ) , : ) ;  %%%  mode  shapes 

end 

end 


-end  of  sol  fv.m 
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F.  eigsn 

% - eigsn. m - 

"o 

function  [V, D] =eigsn (K, M, nmode, tol) 

"o 

%  [V, D]  =  eigsn (K,  M,  nmode) 

%  finds  lowest  eigenvalues  and  eigenvectors  with  mass  normalization 

%  nomde  :  #  of  eigenvalues  to  be  found 

%  defalut  tolerance  =  le-10 

%  [V,  D]  =  eigsn  (K,  M,  nmode,  tol) 

%  tol  :  user  defined  tolerance 

-6 

%  17  September  2007 (jkr) 

"o 

if  (nargin  ==  3) 

tol  =  le-10; 

end 

"o 

%  EIGS (A, K, SIGMA, OPTS)  and  EIGS (A, B, K, SIGMA, OPTS)  specify  options: 

%  OPTS.disp:  diagnostic  information  display  level  [0  |  {1}  |  2] 

%  option^eigs  =  optimset ( ' OPTS . disp ' , 2 ) ; 

"o 

OPTS . disp=0 ; 

[V, D] =eigs (K,M, nmode,  0,  OPTS)  ; 

"6 

%  biorthogonality  condition 

"o 

for  i=l: nmode; 

V ( :  ,  i )  =  V ( : , i) / sqrt (V ( : , i)  '*M*V(:,i)); 

end 

"5 

% - end  of  eigsn. m - 
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G.  MeKe 


MeKe .m 


function  [Me, Ke ] =MeKe (Eid,  DV) 

"o 

%  Generate  element  mass  matrix (Me) and  stiffness  matrix (Ke) 

%  based  on  the  9-node  plate  element  and  corresponding  numbering 

"o 

%  Dof's  are  { U} ' = [phix, phiy , w } ' 


[INPUT] 

-  Eid, Lam, N1 , Elem, Coord_grid, DV, t layer , rho 

-  m_info= [nnel, ndof , edof , ngrid, t_dof , n_elem] ; 


model  information 


%  [OUTPUT] 

%  -  Me:  27x27 

%  -  Ke :  27x27 

"o 

%  [variables] 

%  -  B:  strain  vs.  nodal  displacement  relation  matrix 

%  [ex  ey  rxy] '=z*Bi*[U} ' 

%  -  D:  moment  vs.  curvature  relation  matrix  (cf.  Stress-Strain) 

%  [Mx  My  Mxy] ' = [D] * { kx  ky  ky} ' 

%  (laminate  bending  stiffness  matrix;  Dll , D12 , D16, , , , D66) 

%  -  Qb:  stress  vs.  strain  relation  matrix 

"o 

%  17  September  2007  (jkr) 


global  Coord 
global  s_info 


grid  Elem  Lam  Be; 
m  info; 


%  model  data 

%  specimen  and  model  information 


La=s_info ( 1 ) ;  Lb=s_info (2 ) ;  %  plate  size 

Na=s_info (3) ;  Nb=s_info ( 4 ) ;  %  #  of  elements  in  each  direction,  x  &  y 
Nl=s_info ( 5 ) ;  Ang=s_inf o ( 6 : 5+N1 ) ;  %  #  of  layer  and  lamination  angle 

rho=s_inf o ( 6+N1 ) ;  tlayer=s  info (7+N1) ;  %  density  and  lamina  thickness 


nnel=m_in 
ndof=m_in 
edof=m_in 
ngrid=m  i 
t_dof =m_i 
n  elem=m 


fo  (1)  ; 
fo  (2 )  ; 
fo (3 )  ; 
nfo  ( 4 )  ; 
nfo  ( 5)  ; 
info  (6) , 


number  of  nodes  per  element 
degrees  of  freedom  per  node 
degrees  of  freedom  per  element 
total  number  of  grids  (nodes) 
total  dof 
#  of  elements 


Ke=zeros (edof 
Kb=zeros (edof 
Ks=zeros (edof 
Me=zeros (edof 


,edof);  %  initialization  of  stiffness  matrix 

, edof);  %  init  of  stiffness  matrix  (bending) 

, edof);  %  init  of  stiffness  matrix  (transverse  shear) 

, edof);  %  initialization  of  mass  matrix 


[Db,Ds]=D  matrix (Eid, DV) ;  %  Moment-Curvature  relation  matrix,  D 


Gid=Elem (Eid, 2:10); 
xcoord=Coord_grid (Gid, 1 ) ;  ycoord=Coord_grid (Gid, 2 ) ;  %  grid  coord's 


%  grid  ID's  for  element 


numerical  integration  bending  stiffness  and  mass  (  3x3  integration  ) 


clear  r  wtr  s  wts 
nglxb=3;  nglyb=3; 


%  use  3x3  integration  rule 
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[pt3,  wt3] =glqd2 (nglxb, nglyb) ; 

"6 

for  intx=l: nglxb 
r=pt3 (intx, 1) ; 
wtr=wt3  (intx, 1)  ; 

for  inty=l: nglyb 
s=pt3 (inty, 2) ; 
wts=wt3 (inty, 2) ; 

"o 

%  for  each  integration  points 


%  sampling  pts  &  wts  in  2-D 


%  sampling  point  in  x-axis 
%  weight  in  x-axis 

%  sampling  point  in  y-axis 
%  weight  in  y-axis 


%  compute  shape  functions  and  derivatives  at  integration  point 
[h9, dh9dr , dh9ds ] =shape_iso9 (r , s) ; 

"o 

%  initialize  dh/dx  and  dh/dy 
dh9dx=zeros (l,nnel) ; 
dh9dy=zeros (l,nnel) ; 

"6 

jacob2=zeros (2, 2) ;  %  compute  Jacobian 

for  ii=l:nnel; 

jacob2 (1,1) =jacob2 (1,1) +dh9dr (ii) *xcoord(ii)  ; 
jacob2 (1, 2) =jacob2 (1,2) +dh9dr (ii) *ycoord(ii) ; 
jacob2 (2,1) =jacob2 (2, 1) +dh9ds (ii) *xcoord(ii) ; 
jacob2 (2, 2) =jacob2 (2, 2) +dh9ds (ii) *ycoord(ii) ; 
end 


det j  acob=det ( j  acob2 ) ; 
invjacob=inv ( jacob2) ; 


%  determinant  of  Jacobian 
%  inverse  of  Jacobian  matrix 


for  ii-i:nnel;  %  derivatives  w.r.t  physical  coordinate 

dh9dx (ii) =invjacob (1, 1) *dh9dr (ii) linvjacob (1,2) *dh9ds (ii) ; 
dh9dy (ii) =invjacob (2, 1) *dh9dr (ii) linvjacob (2,2) *dh9ds (ii) ; 

end 

"o 

%  for  bending  stiffness  matrix,  Kb 
-6 

%  Strain-Nodal  displacement  relation  matrix,  Bi 
%  [ex  ey  gxy]'=  z* [Bi ] * [phi_l  phi_2  w] ' 

%  [Bi] = |  dh9dx  0  0 

%  |  . . .  0  dh9dy  0  . . .  | 

%  I  dh9dy  dh9dx  0  I 

"o 

Bi=zeros  (3,27) ; 
for  ii=l:nnel 

Bi  (1, 3* (ii-1) +I)=dh9dx (ii)  ; 

"o 

Bi ( 2 , 3* (ii-1) +2)=dh9dy (ii) ; 

"o 

Bi  (3, 3*  (ii-1) +l)=Bi  (2,3* (ii-1) +2) ; 

Bi (3,  3* (ii-1) +2)=Bi (1,3* (ii-1) +1) ; 

end 

-6 

Kb=Kb+Bi ' *Db*Bi*wtr*wts*det jacob;  %  integration  for  stiffness  matrix 

"o 

%  for  mass  matrix 


%  Displacement  interpolation,  N 
%  [w]  =  [N]  [dof  ]  =  [N]  [  .  .  .  phi_l  phi_2  w  .  .  . 

%  = [ . . .  -zhi  -zhi  hi  ...][.. 

-6 

%  [N]=[-z  -z  1] | hi  0  0  hi  0  0 

%  I  0  hi  0  ...  0  hi  0  . . 

56 


phi_l  phi_2  w  . . . ] ' 

h9  0 
0  h9 


0  I  =  [ z3 ]  [H3 ] 
01 


|0  0  hi 


0  0  hi 


0  0  h9 1 


int (-t/2~t/2)  {  [ z3 ]  '*[z3] } 


(tA3/12) |  1 
I  1 


z33 


I 


|  0  0  12/tA2| 


H3=zeros  (3,27) ; 
for  i=l : nnel 

H3  (1, 3* (i-1) +1) =h9  (i)  ; 

H3 (2,3* (i-1) +2) =h9  (i)  ; 

H3  (3,  3* (i-1) +3) =h9  (i)  ; 

end 

tt=tlayer*Nl; 

z33=tt A3/ 12  * [  1  10; 

1  10; 

0  0  12 / tt A2 ] ; 


%  laminate  thickness 


Me=Me+rho*H3 ' *z33*H3*wtr*wts*det j acob;  %  integration  for  mass  matrix 

"o 

end  %  for  loop  for  integration  in  y 
end  %  for  loop  for  integration  in  x 


Kb= (Kb+Kb' ) / 2; 
Me= (Me+Me ' ) / 2 ; 


%  for  symmetry 
%  for  symmetry 


numerical  integration  trnasverse  shear  stiffness  (  2x2  integration  ) 


clear  r  wtr  s  wts 
nglxs=2;  nglys=2; 

[pt2,wt2] =glqd2 (nglxs ,  nglys ) 


%  use  3x3  integration  rule 
%  sampling  pts  &  wts  in  2-D 


for  intx=l: nglxs 
r=pt2 (intx, 1) ; 
wtr=wt2  (intx, 1)  ; 
for  inty=l: nglys 
s=pt2 (inty, 2) ; 
wts=wt2  (inty, 2) ; 


%  sampling  point  in  x-axis 
%  weight  in  x-axis 

%  sampling  point  in  y-axis 
%  weight  in  y-axis 


%  for  each  integration  points 


%  compute  shape  functions  and  derivatives  at  integration  point 
[h9,  dh9dr , dh9ds ] =shape_iso9  (r , s ) ; 

"o 

%  initialize  dh/dx  and  dh/dy 
dh9dx=zeros (l,nnel) ; 
dh9dy=zeros (l,nnel)  ; 

-O 

jacob2=zeros (2, 2) ;  %  compute  Jacobian 

for  ii=l:nnel; 

jacob2  (1,1) =jacob2 (1,1) +dh9dr (ii) *xcoord(ii) ; 
jacob2 (1,2) =jacob2 (1, 2) +dh9dr (ii) *ycoord(ii)  ; 
jacob2 (2, 1) =jacob2 (2, 1) +dh9ds (ii) *xcoord(ii) ; 
jacob2 (2, 2) =jacob2 (2, 2) +dh9ds (ii) *ycoord(ii)  ; 
end 
"6 

detjacob=det ( jacob2) ;  %  determinant  of  Jacobian 

inv j acob=inv ( j acob2 ) ;  %  inverse  of  Jacobian  matrix 

-6 

for  ii=l:nnel;  %  derivatives  w.r.t  physical  coordinate 

dh9dx (ii) =invjacob (1,1) *dh9dr (ii) linvjacob (1,2) *dh9ds (ii) ; 
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dh9dy (ii) =invjacob (2, 1) *dh9dr (ii) +invjacob (2,2) *dh9ds (ii) ; 


end 

"o 

%  for  shear  stiffness  matrix,  Ks 

"o 


%  Strain-Nodal  displacement  relation  matrix,  Bi 
%  [gxz  gyz]'=  [Bs]*[phi_l  phi_2  w] ' 

%  [ Bs ] = |  ...  -hi  0  dh9dx  7..  | 

%  I  .  .  .  0  -hi  dh9dy  ...  | 

"6 

Bs=zeros  (2,27)  ; 
for  ii=l : nnel 

Bs (1,3* (ii-1) +l)=-h9 (ii) ; 

Bs (1,3* (ii-1) +3)=dh9dx(ii) ; 

-O 

Bs ( 2 , 3* (ii-1) +2)=-h9  (ii) ; 

Bs (2,3* (ii-1) +3)=dh9dy(ii) ; 

end 

"6 

Ks=Ks+Bs ' *Ds*Bs*wtr*wts*det jacob;  %  integration  for  stiffness  matrix 

"o 

end  %  for  loop  for  integration  in  y 
end  %  for  loop  for  integration  in  x 

"o 

Ks= (Ks+Ks ' ) /2 ;  %  for  symmetry 

"o 

-6 

%  total  stiffness  matrix  with  shear  correction  factor, 5/6 


"o 

Ke=Kb+ (5/6) *Ks ; 
-6 
"6 


end  of  MeKe.m 
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H. 


dKedDVi: 


dKe  8Ke 
d(DVt)  dOi 


"o 

"o 

function 

"o 


-  dKedDVi.m  - 

[  dKedDVi, dKedDV2, dKedDV3, dKedDV4 ] =dKedDVi (Eid, DV) 


Generate  derivatives  of  element  stiffness  matrix (Ke)  w.r.t  DVi 
based  on  the  9-node  plate  element  and  corresponding  numbering 

Dof's  are  { U } ' = {phix, phiy , w } ' 


[INPUT] 

-  Eid, Lam, N1 , Elem, Coord_grid, DV, t layer , rho 

-  m_info= [nnel, ndof , edof , ngrid, t  dof , n_elem] ; 


model  information 


[OUTPUT] 


"6 

-  dKedDVi : 

27x27 

derivative 

of 

stiffness 

matirx 

wrt 

DVI 

-6 

-  dKedDV2 : 

27x27 

derivative 

of 

stiffness 

matirx 

wrt 

DV2 

"6 

-  dKedDV3 : 

27x27 

derivative 

of 

stiffness 

matirx 

wrt 

DV3 

"6 

-  dKedDV4 : 

27x27 

derivative 

of 

stiffness 

matirx 

wrt 

DV4 

%  [variables] 

%  -  B:  strain  vs.  nodal  displacement  relation  matrix 

%  [ex  ey  rxy] '=z*Bi*[U} ' 

%  -  D:  moment  vs.  curvature  relation  matrix  (cf.  Stress-Strain) 

%  [Mx  My  Mxy] ' = [D] * { kx  ky  ky} ' 

%  (laminate  bending  stiffness  matrix;  Dll , D12 , D16, , , , D66) 

%  -  Qb:  stress  vs.  strain  relation  matrix 

"o 

%  2  October  2007  (jkr) 


global  Coord_grid  Elem  Lam  Be; 
global  s_info  m_info; 


%  model  data 

%  specimen  and  model  information 


La=s_info ( 1 ) ;  Lb=s_info (2 ) ;  %  plate  size 

Na=s_info (3) ;  Nb=s_info ( 4 ) ;  %  #  of  elements  in  each  direction,  x  &  y 
Nl=s_info ( 5 ) ;  Ang=s_inf o ( 6 : 5+N1 ) ;  %  #  of  layer  and  lamination  angle 

rho=s_inf o ( 6+N1 ) ;  tlayer=s_info ( 7+N1 ) ;  %  density  and  lamina  thickness 


nnel=m_info ( 1 ) ; 
ndof=m_info (2)  ; 
edof=m_info (3)  ; 
ngrid=m_inf o ( 4 )  ; 
t_dof =m_inf o ( 5 )  ; 
n_elem=m_info (6) ; 

"6 

dKedDVl=zeros (edof,  edof)  ; 
dKedDV2=zeros (edof,  edof)  ; 
dKedDV3=zeros (edof,  edof)  ; 
dKedDV4=zeros (edof,  edof)  ; 
"6 

dKbdDVl=zeros (edof,  edof)  ; 
dKbdDV2=zeros (edof,  edof)  ; 
dKbdDV3=zeros (edof,  edof)  ; 
dKbdDV4=zeros (edof,  edof)  ; 

-O 

dKsdDVl=zeros (edof,  edof)  ; 
dKsdDV2=zeros (edof,  edof)  ; 


%  number  of  nodes  per  element 
%  degrees  of  freedom  per  node 
%  degrees  of  freedom  per  element 
%  total  number  of  grids  (nodes) 

%  total  dof 
%  #  of  elements 


"5 

initialization 

of 

matrix 

"5 

initialization 

of 

matrix 

"5 

initialization 

of 

matrix 

"5 

initialization 

of 

matrix 

"5 

initialization 

of 

matrix 

"5 

initialization 

of 

matrix 

"5 

initialization 

of 

matrix 

"o 

initialization 

of 

matrix 

"o 

initialization 

of 

matrix 

"5 

initialization 

of 

matrix 
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dKsdDV3=zeros (edof,edof) ; 
dKsdDV4=zeros (edof , edof ) ; 


%  initialization  of  matrix 
%  initialization  of  matrix 


[dDbdl, dDsdl, dDbd2, dDsd2, dDbd3, dDsd3, dDbd4, dDsd4] =. . . 

D_sen_matrix_wrt_DVi (Eid, DV) ;  %  Derivative  of  D-matrices 

"o 

Gid=Elem (Eid, 2 : 10) ;  %  grid  ID's  for  element 

xcoord=Coord_grid (Gid, 1 ) ;  ycoord=Coord_grid (Gid, 2 ) ;  %  grid  coord's 


numerical  integration  of  derivatives  of  bending  stiffness  (  3x3  int.  ) 


clear  r  wtr  s  wts 
nglxb=3;  nglyb=3; 

[pt3, wt3] =glqd2 (nglxb, nglyb) 

"o 

for  intx=l: nglxb 
r=pt3 (intx, 1) ; 
wtr=wt3 (intx, 1) ; 
for  inty=l: nglyb 
s=pt3 (inty, 2) ; 
wts=wt3 (inty, 2 ) ; 


%  use  3x3  integration  rule 
%  sampling  pts  &  wts  in  2-D 


%  sampling  point  in  x-axis 
%  weight  in  x-axis 

%  sampling  point  in  y-axis 
%  weight  in  y-axis 


%  for  each  integration  points 

"5 

%  compute  shape  functions  and  derivatives  at  integration  point 
[h9, dh9dr , dh9ds ] =shape_iso9 (r , s) ; 

"o 

%  initialize  dh/dx  and  dh/dy 
dh9dx=zeros (l,nnel) ; 
dh9dy=zeros (l,nnel) ; 

"o 

jacob2  =  zeros  (2, 2) ;  %  compute  Jacobian 

for  ii=l:nnel; 

jacob2 (1,1) =jacob2 (1,1) +dh9dr (ii) *xcoord(ii) ; 
jacob2 (1, 2) =jacob2 (1,2) +dh9dr (ii) *ycoord(ii) ; 
jacob2 (2,1) =jacob2 (2, 1) +dh9ds (ii) *xcoord(ii) ; 
jacob2 (2, 2) =jacob2 (2, 2) +dh9ds (ii) *ycoord(ii) ; 
end 


det j  acob=det ( j  acob2 ) ; 
inv j  acob=inv ( j  acob2 ) ; 


determinant  of  Jacobian 
inverse  of  Jacobian  matrix 


for  ii=l:nnel;  %  derivatives  w.r.t  physical  coordinate 

dh9dx (ii) =invjacob (1, 1) *dh9dr (ii) +invjacob (1,2) *dh9ds (ii) ; 
dh9dy (ii) =invjacob (2, 1) *dh9dr (ii) +invjacob (2,2) *dh9ds (ii) ; 

end 

"6 

%  for  bending  stiffness  matrix,  Kb 

"o 

%  Strain-Nodal  displacement  relation  matrix,  Bi 
%  [ex  ey  gxy]'=  z* [Bi ] * [phi_l  phi_2  w] ' 

%  [Bi] = |  dh9dx  0  0  ~| 

%  |  . . .  0  dh9dy  0  . . .  | 

%  I  dh9dy  dh9dx  0  I 

"o 

Bi=zeros  (3,27) ; 
for  ii=l:nnel 

Bi  (1, 3* (ii-1) +1) =dh9dx (ii)  ; 


Bi  ( 2 , 3*  (ii-1) +2)=dh9dy(ii) ; 
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end 


Bi (3,  3* (ii-1) +1) =Bi (2,3* (ii-1) +2)  ; 

Bi (3,  3* (ii-1) +2)=Bi (1,3* (ii-1) +1) ; 


dKbdDVl=dKbdDVl+Bi ' *dDbdl *Bi*wtr*wts*det j  acob ; 
dKbdDV2=dKbdDV2+Bi ' *dDbd2  *Bi*wtr*wts*det j  acob; 
dKbdDV3=dKbdDV3+Bi ' *dDbd3*Bi*wtr*wts*det j  acob ; 
dKbdDV4=dKbdDV4+Bi ' *dDbd4  *Bi*wtr*wts*det j  acob; 
-6 

end  %  for  loop  for  integration  in  y 
end  %  for  loop  for  integration  in  x 
"6 

dKbdDVl= (dKbdDVl+dKbdDVl ' ) / 2 ; 
dKbdDV2= (dKbdDV2+dKbdDV2 ' ) / 2 ; 
dKbdDV3= (dKbdDV3+dKbdDV3 ' ) / 2 ; 
dKbdDV4= (dKbdDV4+dKbdDV4 ' ) / 2 ; 


%  integration 
%  integration 
%  integration 
%  integration 


for  matrix 
for  matrix 
for  matrix 
for  matrix 


%  for  symmetry 
%  for  symmetry 
%  for  symmetry 
%  for  symmetry 


-6 

"o 

%  num  integration  of  derivatives  fo  trans  shear  stiffness  (  2x2  integ  ) 

"6 

"6 

clear  r  wtr  s  wts 


nglxs=2;  nglys=2; 

"6 

use  3x3  integration  rule 

[pt2,wt2] =glqd2 (nglxs , nglys ) ; 

"6 

sampling  pts  &  wts  in  2-D 

for  intx=l: nglxs 

r=pt2 (intx, 1 ) ; 

"6 

sampling  point  in  x-axis 

wtr=wt2 (intx, 1) ; 

-o 

weight  in  x-axis 

for  inty=l: nglys 

s=pt2 (inty, 2) ; 

-6 

sampling  point  in  y-axis 

wts=wt2 (inty, 2) ; 

"6 

weight  in  y-axis 

"6 


%  for  each  integration  points 


%  compute  shape  functions  and  derivatives  at  integration  point 
[h9, dh9dr , dh9ds ] =shape_iso9 (r , s ) ; 

"6 

%  initialize  dh/dx  and  dh/dy 
dh9dx=zeros (l,nnel) ; 
dh9dy=zeros (l,nnel)  ; 

"6 

jacob2=zeros (2, 2) ;  %  compute  Jacobian 

for  ii=l : nnel ; 

jacob2  (1,1) =jacob2 (1,1) +dh9dr (ii) *xcoord(ii) ; 
jacob2 (1, 2) =jacob2 (1, 2) +dh9dr (ii) *ycoord(ii)  ; 
jacob2 (2, 1) =jacob2 (2, 1) +dh9ds (ii) *xcoord(ii) ; 
jacob2 (2, 2) =jacob2 (2, 2) +dh9ds (ii) *ycoord(ii) ; 
end 
"6 

detjacob=det ( jacob2) ;  %  determinant  of  Jacobian 

inv j acob=inv ( j acob2 ) ;  %  inverse  of  Jacobian  matrix 

"o 

for  ii=l:nnel;  %  derivatives  w.r.t  physical  coordinate 

dh9dx (ii) =invjacob (1,1) *dh9dr (ii) linvjacob (1,2) *dh9ds (ii) ; 
dh9dy (ii) =invjacob (2,1) *dh9dr (ii) linvjacob (2,2) *dh9ds (ii) ; 

end 

"o 

%  for  shear  stiffness  matrix,  Ks 
-6 


%  Strain-Nodal  displacement  relation  matrix,  Bi 
%  [gxz  gyz]'=  [Bs]*[phi_l  phi_2  w] ' 
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%  [Bs]=|  ...  -hi  0  dh9dx  ... 

%  I  . . .  0  -hi  dh9dy 

"o 

Bs=zeros  (2,27)  ; 
for  ii=l : nnel 

Bs (1, 3* ( i i  —  1 ) +1 ) =-h9 ( ii ) ; 

Bs (1,3* ( i i — 1 ) +3 ) =dh9dx ( ii ) ; 

-6 

Bs ( 2 , 3* (ii-1) +2)=-h9  (ii) ; 

Bs (2,3* (ii-1) +3)=dh9dy(ii) ; 

end 


dKsdDVl=dKsdDVl+Bs ' *dDsdl *Bs*wtr*wts*det j  acob ; 
dKsdDV2=dKsdDV2+Bs ' *dDsd2  *Bs*wtr*wts*det j  acob; 
dKsdDV3=dKsdDV3+Bs ' *dDsd3*Bs*wtr*wts*det j  acob; 
dKsdDV4=dKsdDV4+Bs ' *dDsd4  *Bs*wtr*wts*det j  acob ; 

"o 

end  %  for  loop  for  integration  in  y 
end  %  for  loop  for  integration  in  x 

"o 

dKsdDVl= (dKsdDVl+dKsdDVl ' ) / 2 ; 
dKsdDV2= (dKsdDV2+dKsdDV2 ' ) / 2 ; 
dKsdDV3= (dKsdDV3+dKsdDV3 ' ) / 2 ; 
dKsdDV4= (dKsdDV4+dKsdDV4 ' ) / 2 ; 


integration 

integration 

integration 

integration 


for  matrix 
for  matrix 
for  matrix 
for  matrix 


%  for  symmetry 
%  for  symmetry 
%  for  symmetry 
%  for  symmetry 


%  total  stiffness  matrix  with  shear  correction  factor, 5/6 


dKedDVl=dKbdDVl+ (5/6) *dKsdDVl; 
dKedDV2=dKbdDV2+ (5/6) *dKsdDV2; 
dKedDV3=dKbdDV3+ (5/6) *dKsdDV3; 
dKedDV4=dKbdDV4+ (5/6) *dKsdDV4; 

"o 

% - end  of  dKedDVi.m 
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I.  D  matrix 


% - D_matrix.m 

"o 

function  [Db, Ds ] =D_matrix (Eid, DV) 


"o 


"o 

D  matrices  of 

each 

element 

"o 

"o 

[INPUT] 

"6 

-  Eid, 

DV 

"o 

"o 

[OUTPUT] 

-o 

-  Db  = 

|  Dll 

D12 

D16  | 

"bending  stiffness" 

"o 

|  D12 

D22 

D2  6  I 

"o 

|  D16 

D2  6 

D66  | 

"o 

-  Ds  = 

|  D55 

D54 

1 

"transverse  shear  stiffness" 

"o 

|  D54 

D44 

1 

"6 

%  17  September  2007  (jkr) 

"o 

global  Coord_grid  Elem  Lam  Be;  %  model  data 

global  s_info  m_info;  %  specimen  and  model  information 

"o 

Nl=s_info  (5 )  ; 

"o 

Lam_=Lam(Nl* (Eid-1) +1 ;Nl*Eid, :,:,:); 
layer_id=Lam_ (1 :N1,  2)  ; 
layer_t=Lam_ ( 1 : N1 , 3 )  ; 
layer_angle=Lam_ ( 1 : N1 , 4 ) ; 

"o 

E1=DV ( 1 ) ; 

E2=DV (2 ) ; 
vl2=DV (3) ; 

G12=DV ( 4 ) ; 
v21=vl2*E2/El; 

"O 

tt=sum (layer_t) ;  %  laminate  thickness 

"o 

%  in-plane  property 

"o 

Q11=E1/ (1-v12*v21 ) ; 

Q22=E2/ (1-v12*v21 ) ; 

Q12=vl2*E2/ (1-v12*v21) ; 

Q  6  6  =G 1 2 ; 

"O 

%  trnasverse  shear  properties  are  expressed  in  Q66  and  0 . 5* (Q11-Q12) 
%  (assumed  transversely  isotropic) 

-6 

z  ( 1 ) =-tt/ 2 ; 
zcoord=z ( 1 ) ; 
for  ii=l:Nl 

zcoord=zcoord+layer_t (ii)  ; 
z (ii+1) =zcoord; 

end 

-6 

Db=zeros (3,3)  ; 

Ds=zeros  (2,2)  ; 

%B=zeros (3, 3)  ; 

%A=zeros  (3,3)  ; 
for  ii=l:Nl 
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o\o  o\o 


theta=layer_angle (ii) *pi/180; 
cs=cos (theta) ;  sn=sin (theta) ; 

Qbll=Qll*csA4+2* (Q12+2*Q66) *snA2*cs A2+Q22 *snA 4 ; 

Qbl2=  (Q11+Q22-4*Q66) *snA2*csA2+Q12* (snA4  +  csA4) ; 
Qb22=Qll*snA4+2* (Q12+2*Q66) *snA2*cs A2+Q22 *csA 4 ; 

Qbl 6=  (Q11-Q12-2*Q66) *sn*csA3+  (Q12-Q22+2*Q66) *snA3*cs 
Qb2  6= (Q11-Q12-2*Q66) *snA3*cs+ (Q12-Q22+2*Q66) *sn*csA3 
Qb6 6= (Q11+Q22-2*Q12-2*Q66) * snA2 *cs A2+Q6 6*  (snA4  +  csA4) 

Qs55=Q6  6*cs A2  +  0 .5* (Q11-Q12) *snA2; 
Qs54=-Q66*cs*sn+0.5* (Q11-Q12) *cs*sn; 

Qs44=Q66*snA2+0 .5* (Q11-Q12) *csA2; 

Qb= [Qbl 1  Qbl 2  Qbl6; 

Qbl2  Qb22  Qb2 6 ; 

Qbl 6  Qb2 6  Qb66] ; 

Qs=[Qs55  Qs54; 

Qs54  Qs44 ] ; 

Db=Db  +  l/3*Qb* (z (ii+1) A3-z  (ii) A3) ; 

Ds=Ds+Qs* (z (ii  +  1) -z  (ii)  ) ; 

B=B+l/2*Qb* (z (ii  +  1) A2-z  (ii) A2)  ; 

A=A+Qb* ( z (ii  +  1 ) -z  (  ii ) )  ; 

end 

% - end  of  D  matrix. m - 
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J.  D  sen  matrix  wrt  DVi 


% -  D_sen_matr ix_wrt_DVi .m  - 

"6 

function  [dDbdl , dDsdl , dDbd2 , dDsd2 , dDbd3 , dDsd3 , dDbd4 , dDsd4 ] = . . . 
D_sen_matrix_wrt_DVi (Eid, DV) 

"6 

%  Sensitivity  of  D_matrices (Db  and  Ds)  w.r.t.  DV1=E1 , DV2=E2 , DV3=vl2 , DV3=G12 


[INPUT] 


"6 

-  Eid,  N1 

,  Lam,  DV 

-6 

[OUTPUT] 

"6 

-  dDbdi  = 

d  1 

Dll 

D12  D16 

"6 

-  1 

D12 

D22  D26 

"6 

d (DVi )  | 

D16 

D26  D66 

% 

"6 

-  dDsdi  = 

d  I 

D55 

D54  | 

"O 

"6 

1 

d  (DVi )  | 

D54 

1 

D44  | 

"o 

%  l 

October  2007 

( j  kr) 

global  Coord_grid  Elem  Lam  Be;  %  model  data 

global  s_info  m_info 
"6 

La=s_info ( 1 ) ;  Lb=s_info (2 ) ;  %  plate  size 

Na=s_info (3) ;  Nb=s_info ( 4 ) ;  %  #  of  elements  in  each  direction,  x  &  y 
Nl=s_info ( 5 ) ;  Ang=s_inf o ( 6 : 5+N1 ) ;  %  #  of  layer  and  lamination  angle 

rho=s_inf o ( 6+N1 ) ;  tlayer=s_info ( 7+N1 ) ;  %  density  and  lamina  thickness 


nnel=m_info ( 1 ) ; 
ndof=m_info (2)  ; 
edof=m_info (3)  ; 
ngrid=m_inf o ( 4 )  ; 
t_dof =m_inf o ( 5 )  ; 
n  elem=m  info  (6); 


%  number  of  nodes  per  element 
%  degrees  of  freedom  per  node 
%  degrees  of  freedom  per  element 
%  total  number  of  grids  (nodes) 

%  total  dof 
%  #  of  elements 


Lam_=Lam(Nl* (Eid-1) +1 ;Nl*Eid, 
layer_id=Lam_ (1 :N1, 2) ; 
layer_t=Lam_ ( 1 : N1 , 3 )  ; 
layer_angle=Lam_ ( 1 : N1 , 4 ) ; 


E1=DV(1) ; 

E2=DV (2 ) ; 
vl2=DV (3) ; 

G12=DV ( 4 ) ; 
v2 1  =v  1 2  *E  2  /E 1 ; 

"6 

tt=sum (layer_t) ;  %  laminate  thickness 

"6 

%  on-axis  properties 
"6 

%  in-plane 

%  Q11=E1/ (1-v12*v21) ; 

%  Q22=E2/ (1-v12*v21)  ; 

%  Q12=vl2*E2/ (1-v12*v21)  ; 

%  Q  6  6  =G 1 2 ; 

"6 

%  trnasverse  shear  (assumed  transversely  isotropic) 
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Q55=Q66 

Q44=0 .5* (Q11-Q12) 


"O 
"6 

"o 

%  Derivatives  wrt  DV1;  d (Qll) /d (El) ,  d (Q22 ) /d (El )  ,  ,  , 

"o 

dQlldl=l/ (1-v12*v21) -vl2*v21/ (1-v12*v21) A2; 
dQ22dl=-v21A2/ (1-v12*v21) A2; 
dQ12dl=-vl2*v21A2/ (1-v12*v21) A2; 
dQ6  6dl  =  0 ; 

"o 

%  Derivatives  wrt  DV2 ;  d (Qll ) /d (E2 ) ,  d (Q22) /d (E2 ) , , , 

"o 

dQlld2=vl2A2/  (1-v12*v21) A2; 
dQ22d2=l/ (1-v12*v21)+v12*v21/ (1-v12*v21) A2; 
dQ12d2=vl2/ (1-v12*v21 ) +vl2A2*v21/ (1-v12*v21) A2; 
dQ6  6d2  =  0 ; 

"6 

%  Derivatives  wrt  DV3;  d (Qll ) /d (vl2 ) ,  d (Q22 ) /d (vl2) , , , 

"6 

dQlld3=2*vl2*E2/ (1-v12*v21) A2; 
dQ22d3=2*v21*E2/ (1-v12*v21) A2; 

dQ12d3=E2/ ( 1-v12*v2 1 ) +2*vl2*v21*E2/  (1-v12*v21) A2; 
dQ6  6d3=0 ; 

"6 

%  Derivatives  wrt  DV1;  d (Qll ) /d (E2 ) ,  d (Q22 ) /d (E2 ) , , , 

"6 

dQl ld4=0 ; 
dQ22d4=0; 
dQ12d4=0; 
dQ6  6d4  =  l ; 

"6 

z  ( 1 ) =-tt/ 2 ; 
zcoord=z ( 1 ) ; 
for  ii=l:Nl 

zcoord=zcoord+layer_t (ii)  ; 
z (ii+1) =zcoord; 

end 

-6 

dDbdl=zeros  (3, 3) ;  dDsdl=zeros  (2 , 2 ) ;  %  initialize 

dDbd2=zeros (3, 3) ;  dDsd2=zeros (2 , 2 ) ;  %  initialize 

dDbd3=zeros (3, 3) ;  dDsd3=zeros (2 , 2 ) ;  %  initialize 

dDbd4=zeros (3, 3) ;  dDsd4=zeros (2 , 2 ) ;  %  initialize 

-6 

%  sensitivity  wrt  DV1 

"6 

for  ii=l:Nl 

theta=layer_angle (ii) *pi/180; 
cs=cos (theta) ;  sn=sin (theta) ; 

"o 

dQblldl=dQlldl*csA4+2* (dQ12dl+2*dQ66dl ) * snA2 *csA2+dQ22dl *snA 4 ; 

dQbl2dl= (dQlldl+dQ22dl-4*dQ66dl) *snA2*cs A2+dQ12dl* (snA4+csA4) ; 

dQb22dl=dQlldl*snA4+2* (dQ12dl+2*dQ66dl ) * snA2 *cs A2+dQ22dl *cs A4 ; 

dQbl 6dl= (dQlldl-dQ12dl-2*dQ66dl) *sn*csA3+ (dQ12dl-dQ22dl+2*dQ66dl ) *snA3*cs; 

dQb2  6dl= (dQlldl-dQ12dl-2*dQ6  6dl) *snA3*cs+ (dQ12dl-dQ22dl  +  2*dQ66dl ) *sn*csA3; 

dQb6  6dl= (dQlldl+dQ22dl-2*dQ12dl-2*dQ66dl) *snA2*csA2+dQ66dl* (snA4+csA4) ; 

"6 

dQs55dl=dQ66dl*csA2+0 . 5* (dQl ldl-dQ12dl ) *snA2; 
dQs54dl=-dQ6  6dl*cs*sn+0 . 5* (dQl ldl-dQ12dl ) *cs*sn; 
dQs44dl=dQ66dl*snA2+0 . 5* (dQl ldl-dQ12dl ) *csA2; 

"o 

dQbdl= [dQblldl  dQbl2dl  dQbl6dl; 
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dQbl2dl  dQb22dl  dQb26dl; 
dQbl 6dl  dQb2  6dl  dQb66dl] ; 


dQsdl= [dQs55dl  dQs54dl; 

dQs54dl  dQs44dl] ; 

"6 

dDbdl=dDbdl  +  l/3*dQbdl* (z (ii+1)  A3-z (ii)  A3)  ; 
dDsdl=dDsdl+dQsdl *(z(ii+l)-z(ii)); 
end 


%  sensitivity  wrt  DV2 


for  ii=l:Nl 

theta=layer_angle (ii) *pi/180; 
cs=cos (theta) ;  sn=sin (theta) ; 

"6 

dQbl ld2=dQlld2*csA 4+2* (dQ12d2+2*dQ66d2 ) * snA2 *cs A2+dQ22d2 *snA 4 ; 

dQbl2d2= (dQl ld2+dQ22d2-4  *dQ6  6d2 ) *snA2*csA2+dQ12d2* (snA4  +  csA4) ; 

dQb22d2=dQlld2*snA4  +  2* (dQ12d2  +  2*dQ66d2 ) * snA2 *cs A2+dQ22d2 *csA 4 ; 

dQbl 6d2= (dQl ld2-dQ12d2-2  *dQ6  6d2 ) *sn*csA3+ (dQ12d2-dQ22d2  +  2*dQ66d2 ) *snA3*cs; 

dQb2  6d2= (dQlld2-dQ12d2-2*dQ66d2) *snA3*cs+ (dQ12d2-dQ22d2  +  2*dQ66d2 ) *sn*csA3; 

dQb6  6d2= (dQl ld2+dQ22d2-2 *dQ12d2-2*dQ66d2 ) *snA2*csA2+dQ66d2* (snA4+csA4) ; 

"o 

dQs55d2=dQ66d2*csA2+0 . 5* (dQl ld2-dQ12d2 ) *snA2; 
dQs54d2=-dQ6  6d2*cs*sn+0 . 5* (dQl ld2-dQ12d2 ) *cs*sn; 
dQs44d2=dQ66d2*snA2+0 . 5* (dQl ld2-dQ12d2 ) *csA2; 

"6 

dQbd2  = [ dQb 1 1 d2  dQbl2d2  dQbl6d2; 

dQbl2d2  dQb22d2  dQb26d2; 
dQbl 6d2  dQb2 6d2  dQb66d2]; 

"o 

dQsd2= [dQs55d2  dQs54d2; 

dQs54d2  dQs44d2]; 

-6 

dDbd2=dDbd2  +  l/3*dQbd2* (z (ii  +  1) A3-z (ii) A3)  ; 
dDsd2=dDsd2+dQsd2* (z (ii+1) -z (ii) ) ; 
end 


%  sensitivity  wrt  DV3 


for  ii=l:Nl 

theta=layer_angle (ii) *pi/180; 
cs=cos (theta) ;  sn=sin (theta) ; 

"o 

dQbl ld3=dQlld3*csA 4  +  2* (dQ12d3  +  2*dQ66d3) * snA2 *cs A2+dQ22d3*snA 4  ; 

dQbl2d3= (dQl ld3+dQ22d3-4 *dQ6 6d3) *snA2*csA2+dQ12d3* (snA4  +  csA4)  ; 

dQb22d3=dQlld3*snA4+2* (dQ12d3+2*dQ66d3) * snA2 *cs A2+dQ22d3*cs A4 ; 

dQbl 6d3= (dQl ld3-dQ12d3-2 *dQ6 6d3) *sn*csA3+ (dQ12d3-dQ22d3  +  2*dQ66d3)  *snA3*cs; 

dQb2  6d3= (dQlld3-dQ12d3-2*dQ66d3) *snA3*cs+ (dQ12d3-dQ22d3  +  2*dQ66d3) *sn*csA3; 

dQb6  6d3= (dQlld3+dQ22d3-2*dQ12d3-2*dQ66d3) *snA2*csA2+dQ66d3* (snA4+csA4) ; 

"o 

dQs55d3=dQ66d3*csA2+0 . 5* (dQl ld3-dQ12d3 ) *snA2; 
dQs54d3=-dQ66d3*cs*sn+0.5*(dQlld3-dQ12d3)*cs*sn; 
dQs4  4d3=dQ66d3*snA2  +  0 . 5* (dQl ld3-dQ12d3 )  *csA2; 

"6 

dQbd3= [ dQb 1 1 d3  dQbl2d3  dQbl6d3; 

dQbl2d3  dQb22d3  dQb26d3; 
dQbl 6d3  dQb2 6d3  dQb66d3]; 

"o 

dQsd3= [dQs55d3  dQs54d3; 
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dQs54d3  dQs44d3]; 


dDbd3=dDbd3  +  l/3*dQbd3* (z (ii+1) A3-z (ii) A3)  ; 
dDsd3=dDsd3+dQsd3* (z (ii  +  1) -z  (ii) ) ; 
end 

"6 

"6 

%  sensitivity  wrt  DV4 


for  ii=l : N1 

theta=layer_angle (ii) *pi/180; 
cs=cos (theta) ;  sn=sin (theta) ; 

"6 

dQblld4=dQlld4*csA4+2* (dQ12d4+2*dQ66d4 ) * snA2 *cs A2+dQ22d4 *snA 4 ; 

dQbl2d4= (dQlld4+dQ22d4-4*dQ6  6d4) *snA2*csA2+dQ12d4* (snA4  +  csA4)  ; 

dQb22d4=dQlld4*snA4+2* (dQ12d4+2*dQ66d4 ) * snA2 *cs A2+dQ22d4 *csA 4 ; 

dQbl 6d4= (dQlld4-dQ12d4-2*dQ66d4) *sn*csA3+ (dQ12d4-dQ22d4+2*dQ66d4 ) *snA3*cs; 

dQb2  6d4= (dQl ld4-dQ12d4-2 *dQ6 6d4 ) *snA3*cs+ (dQ12d4-dQ22d4  +  2*dQ66d4 ) *sn*csA3; 

dQb6  6d4= (dQl ld4+dQ22d4-2 *dQ12d4-2*dQ66d4 ) *snA2*csA2+dQ66d4* (snA4+csA4) ; 

"o 

dQs55d4=dQ66d4*csA2+0 . 5* (dQl ld4-dQ12d4 ) *snA2; 
dQs54d4=-dQ6  6d4*cs*sn+0 . 5* (dQl ld4-dQ12d4 ) *cs*sn; 
dQs44d4=dQ66d4*snA2+0. 5* (dQl ld4-dQ12d4 ) *csA2; 

"o 

dQbd4= [dQblld4  dQbl2d4  dQbl6d4; 

dQbl2d4  dQb22d4  dQb26d4; 
dQbl 6d4  dQb2 6d4  dQb66d4]; 

"6 

dQsd4= [dQs55d4  dQs54d4; 

dQs54d4  dQs44d4]; 

-6 

dDbd4=dDbd4  +  l/3*dQbd4* (z (ii  +  1) A3-z (ii) A3)  ; 
dDsd4=dDsd4+dQsd4* (z (ii  +  1) -z (ii)  )  ; 
end 

% - end  of  D  sen  matrix  wrt  DVi.m - 
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K.  seneigenvalue 


% -  sen_eigvalue . m  - 

"o 

function  dn_lamdn_DV=sen_eigvalue (n__DV) 

"o 

'o 

%  sensitivity  of  eigenvalues  w.r.t  design  variables 

"O 

-6 

global  DV_r; 

global  Coord_grid  Elem  Lam  Be;  %  model  data 

global  s_info  m_info;  %  specimen  and  model  information 

global  nmode; 

global  f  exp  v  exp  mac;  %  experimental  results,  weighting 

global  f_ana  v_ana  mshp_gid  lambda  Phi;  %  defined  in  the  routine,  pi_fv_r2 


La=s_info ( 1 ) ;  Lb=s_info  (2 ) ;  %  plate  size 

Na=s_info (3 ) ;  Nb=s_info ( 4 ) ;  %  #  of  elements  in  each  direction,  x  &  y 
Nl=s_info ( 5 ) ;  Ang=s_inf o ( 6 : 5+N1 ) ;  %  #  of  layer  and  lamination  angle 

rho=s_inf o ( 6+N1 ) ;  tlayer=s_info ( 7+N1 ) ;  %  density  and  lamina  thickness 


nnel=m_info ( 1 ) ; 
ndof=m_info (2)  ; 
edof=m_info (3)  ; 
ngrid=m_inf o ( 4 )  ; 
t_dof=m_info ( 5)  ; 
n_elem=m_info (6) ; 

"o 

DV=n  DV.*DV  r; 


%  number  of  nodes  per  element 
%  degrees  of  freedom  per  node 
%  degrees  of  freedom  per  element 
%  total  number  of  grids  (nodes) 

%  total  dof 
%  #  of  elements 


%  assemble  of  stiffness  sensitivity  matrix 


"o 


dKdDVl  =  zeros  (t 

dof. 

t 

dof)  ; 

"6 

initialization 

of 

stiffness 

derivative 

wrt 

DV1 

dKdDV2  =  zeros  (t 

dof. 

t 

dof)  ; 

"6 

initialization 

of 

stiffness 

derivative 

wrt 

DV2 

dKdDV3=zeros  (t 

dof. 

t 

dof)  ; 

-6 

initialization 

of 

stiffness 

derivative 

wrt 

DV3 

dKdDV4  =  zeros  (t 

dof. 

t 

dof)  ; 

-6 

initialization 

of 

stiffness 

derivative 

wrt 

DV4 

-6 

for  Eid=l:n  elem; 


-5 

[dKedDVl, dKedDV2, dKedDV3, dKedDV4] =dKedDVi (Eid, DV) ; 


gids=Elem (Eid, 2 : 10) ;  %  grid  id's  of  the  element 

idx=feeldof (gids, nnel, ndof ) ; %  system  dof's  of  the  element 


dKdDVl=feasmbll (dKdDVl , dKedDVl , idx) ; 
dKdDV2=feasmbll (dKdDV2, dKedDV2, idx) ; 
dKdDV3=feasmbll (dKdDV3, dKedDV3, idx) ; 
dKdDV4=feasmbll (dKdDV4 , dKedDV4 , idx) ; 


"5 

assemble 

of 

matrix 

for 

DV1 

-6 

assemble 

of 

matrix 

for 

DV2 

"o 

assemble 

of 

matrix 

for 

DV3 

"o 

assemble 

of 

matrix 

for 

DV4 

end 

"6 

%  apply  boundary  condition  (cantilever  plate  -  LHS  clamped) 


[nr , nc] =size (Be) ; 
nbc=nr* (nc-1 )  ;  % 

dKadDVl=dKdDVl (nbc+1 ; t 
dKadDV2=dKdDV2 (nbc+1 : t 
dKadDV3=dKdDV3 (nbc+1 : t 
dKadDV4=dKdDV4 (nbc+1 : t 


total  #  of  dof  constrained 
_dof, nbc+1 ; t_dof) ;  %  matrix 

_dof, nbc+1 ; t_dof) ;  %  matrix 

_dof, nbc+1 : t_dof) ;  %  matrix 

dof, nbc+1 :t  dof);  %  matrix 


patition  (apply 
patition (apply 
patition  (apply 
patition  (apply 


BC'S) 

BC'S) 

BC'S) 

BC'S) 
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%  normalized  sensitivity  of  eigenvalue  for  each  mode 
-6 

dn_lamdn_DV=zeros (length (n_DV) , nmode) ; 

"o 

for  m=l: nmode 

dn_lamdn_DV (1 , m) =Phi ( : , m) ' *dKadDVl*Phi ( : , m) *DV_r ( 1 
dn_lamdn_DV (2 , m) =Phi ( : , m) ' *dKadDV2*Phi ( : , m) *DV_r ( 2 
dn_lamdn_DV (3 , m) =Phi ( : , m) ' *dKadDV3*Phi ( : , m) *DV_r ( 3 
dn_lamdn_DV (4,m)=Phi(:,m) ' *dKadDV4*Phi ( : , m) *DV_r ( 4 

end 

"o 

% - end  of  sen_eigvalue  .m - 


) *l/lambda (m) 
) *l/lambda (m) 
) *l/lambda (m) 
) *l/lambda (m) 
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L.  shape_iso9 


% -  shape_iso9.m - 

"o 

function  [h9, dh9dr, dh9ds] =shape_iso9 (r, s) 

"o 

%  Returns  the  value  of  following  functions  at  corresponding  node (r, s) 


'o 

h9 : 

shape  function  of  9-node  element 

'o 

dh9dr : 

d  ( h 9 ) /dr 

-6 

dh9ds : 

d  ( h 9 ) /ds 

~Q 

%  node  numbering  as  shown  below: 


-6 


"6 

A  s 

~o 

1 

'o 

o- 

— 

- o - 

- o 

'o 

l 

4 

7 

1 

3 

“6 

l 

1 

'o 

o 

o 

o 

- >  r 

'o 

1 

8 

9 

1 

6 

“6 

1 

1 

'o 

o- 

— 

- o - 

- o 

'o 

1 

5 

2 

"o 

%  Refer  to  the  Bathe  or  Kwon  and  Bang 
%  No:  node  number  (1  to  9) 

%  r,s:  node  coordinates 

"o 

%  15  July  2007  (jkr) 

"o 

h9=zeros (1, 9) ;  dh9dr=zeros ( 1 , 9) ;  dh9ds=zeros ( 1 ,  9 )  ; 

"o 

%  Shape  function 
-6 

h9 (1) =0 .25* (r*r-r) * (s*s-s)  ; 
h9 (2) =0 .25* (r*r+r ) * (s*s-s)  ; 
h9 (3)  =0 .2  5* (r*r+r ) *  (s*s  +  s) ; 
h9 (4)  =0 .25* (r*r-r) * (s*s+s) ; 
h9 (5)  =-0. 5* (r*r-l) * (s*s-s)  ; 
h9 (6)  =-0. 5* (r*r+r ) *  (s*s-l) ; 
h9(7)=-0.5* (r*r-l ) * (s*s+s) ; 
h9 (8)  =-0. 5* (r*r-r ) *  (s*s-l) ; 
h9 ( 9)  —  (r*r-l) * (s*s-l)  ; 

"6 

%  Derivative  of  shape  function  w.r.t  "r" 

"6 

dh9dr ( 1 ) =0 . 25* ( 2*r-l ) * ( s* s-s) ; 
dh9dr (2 ) =0 . 25* ( 2*r+l ) * ( s* s-s)  ; 
dh9dr (3 ) =0 . 25* ( 2*r+l ) * ( s* s  +  s)  ; 
dh9dr ( 4 ) =0 . 25* ( 2*r-l ) *  ( s* s  +  s) ; 
dh9dr (5) =-0.5* ( 2*r ) * ( s* s-s)  ; 
dh9dr ( 6 ) =-0 . 5* ( 2*r+l ) *  ( s* s-1 )  ; 
dh9dr (7)=-0.5*(2*r)*(s*s+s); 
dh9dr (8) =-0 . 5* (2*r-l) * (s*s-l)  ; 
dh9dr (9)  =  (2*r) * (s*s-l)  ; 

'O 

%  Derivative  of  shape  function  w.r.t  "s" 

'O 

dh9ds ( 1 ) =0 . 25* (r*r-r) *  (2*s-l)  ; 
dh9ds (2 ) =0 . 25* (r*r+r) *  (2*s-l)  ; 
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dh9ds (3) =0 . 25* (r*r+r) * (2*s  +  l)  ; 
dh9ds ( 4 ) =0 . 25* (r*r-r) * (2*s+l) ; 
dh9ds (5) =-0 . 5* (r*r-l) *  ( 2 * s  —  1 ) ; 
dh9ds (6) =-0 . 5* (r*r+r) *  (2*s) ; 
dh9ds (7 ) =-0 . 5* ( r*r-l) * ( 2*  s  +  1 ) ; 
dh9ds (8) =-0 . 5* (r*r-r) *  (2*s)  ; 
dh9ds ( 9)  =  (r*r-l ) *  (2*s)  ; 

"6 

% - end  of  shape_iso9.m 
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M 


feasmbll 


feasmbll .m 


function  [kk] =feasmbll (kk, k, index) 


Purpose : 

Assembly  of  element  matrices  into  the  system  matrix 
Synopsis : 

[ kk] = feasmbll ( kk, k,  index) 

Variable  Description: 
kk  -  system  matrix 
k  -  element  matri 

index  -  d.o.f.  vector  associated  with  an  element 


edof  =  length  (index) ; 
for  i=l : edof 


ii=index ( i ) ; 
for  j  =1 : edof 

j j=index(j)  ; 

kk (ii, j  j ) =kk (ii, j  j ) +k (i, j ) ; 

end 


end 

-6 

% - end  of  feasmbll. m 
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N, 


feeldof 


feeldof .m 


function  index=feeldof (nd, nnel,  ndof ) 


Purpose : 

Compute  system  dofs  associated  with  each  element 
Synopsis : 

[ index] =fee ldo f (nd,  nnel, ndof) 

Variable  Description: 

index  -  system  dof  vector  associated  with  element  "iel" 
iel  -  element  number  whose  system  dofs  are  to  be  determined 
nd  -  grid  numbers  of  the  element 
nnel  -  number  of  nodes  per  element 
ndof  -  number  of  dofs  per  node 


k=0  ; 

for  i=l : nnel 

start  =  (nd (i) -1) *ndof  ; 
for  j  =1 : ndof 
k=k+l ; 

index (k) =start+j ; 

end 


end 

~o 

% - end  of  feeldof. m 
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o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\°  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  hh  o\o  o\o 


O.  feglqdl 


feglqdl .m 


unction 


[point 1 , weightl ] =f eglqdl (ngl ) 


Purpose : 

determine  the  integration  points  and  weighting  coefficients 
of  Gauss-Legendre  quadrature  for  one-dimensional  integration 


Synopsis : 

[pointl,weightl]=feglqdl (ngl) 


Variable  Description: 

ngl  -  number  of  integration  points 

pointl  -  vector  containing  integration  points 

weightl  -  vector  containing  weighting  coefficients 


initialization 

pointl=zeros (ngl , 1 ) ; 
weightl  =  zeros (ngl,  1)  ; 

find  corresponding  integration  points  and  weights 

if  ngl==l  %  1-point  quadrature  rule 

pointl ( 1 ) =0 . 0 ; 
weightl  ( 1 ) =2 . 0 ; 

"o 

elseif  ngl==2  %  2-point  quadrature  rule 

pointl (1) =-l/sqrt (3) ;  %  pointl ( 1 ) =-0 . 57 73502 6 91 8 9 62 6 
pointl (2) tt-pointl (1 ) ; 
weightl ( 1 ) =1 . 0 ; 
weightl  (2 ) =weightl ( 1 ) ; 

"o 

elseif  ngl==3  %  3-point  quadrature  rule 

pointl (1) =-sqrt (3/5) ;  %pointl (1 ) =-0 . 774596669241483 
pointl (2) =0 . 0; 
pointl (3) =-pointl  (1); 

weightl  (1) =5/9;  %  weightl  ( 1 ) =0 . 55555555555555 6 
weightl (2 ) =8/9;  %  weightl  ( 1 ) =0 . 88 88 88 88 88 88 88 9 
weightl (3 ) =weightl ( 1) ; 

"o 

elseif  ngl==4  %  4-point  quadrature  rule 

pointl (1) =-0. 8 611 3631 15 94 053; 
pointl (2) =-0. 339981043584856; 
pointl (3) =-pointl  (2) ; 
pointl (4) =-pointl  (1)  ; 
weightl (1 ) =0 . 347854845137454 ; 
weightl  (2)=0.6521451548  62  546; 
weightl (3 ) =weightl (2 ) ; 
weightl ( 4 ) =weightl ( 1)  ; 

"o 

else  %  5-point  quadrature  rule 

pointl (1) =-0. 90 61 7 984 5938664; 
pointl (2) =-0. 5384 69310105683; 
pointl (3) =0 . 0; 
pointl ( 4 ) =-pointl  (2 ) ; 
pointl (5) =-pointl  (1) ; 
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o\°  o\o 


weightl ( 1 ) =0 . 23 692 688 505618 9 ; 
weightl (2 ) =0 . 47 8 628 67 04 9936 6 ; 
weightl  ( 3 ) =0 . 5688 88 88 88 88 88 9 ; 
weightl ( 4 ) =weightl  ( 2 ) ; 
weightl ( 5 ) =weightl  ( 1 ) ; 

end 

end  of  feglqdl.m 
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o\°  o\o  o\o  o\o  o\o  o\o  o\°  o\o  o\o  o\°  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  hh  o\o  o\o 


P.  feglqd2 


-  feglqd2 .m 

ion  [point2, weight 2 ] =f eglqd2 (nglx,  ngly ) 


Purpose : 

determine  the  integration  points  and  weighting  coefficients 
of  Gauss-Legendre  quadrature  for  two-dimensional  integration 

Synopsis : 

[point 2 , weight 2 ] =feglqd2 (nglx, ngly) 

Variable  Description: 

nglx  -  number  of  integration  points  in  the  x-axis 
ngly  -  number  of  integration  points  in  the  y-axis 
point2  -  vector  containing  integration  points 
weight2  -  vector  containing  weighting  coefficients 


determine  the  largest  one  between  nglx  and  ngly 

if  nglx  >  ngly 
ngl=nglx; 
else 

ngl=ngly; 

end 

initialization 

point2=zeros (ngl ,2) ; 
weight2  =  zeros (ngl,  2)  ; 

find  corresponding  integration  points  and  weights 

[pointx,  weightx] =f eglqdl (nglx)  ;  %  quadrature  rule  for  x-axis 

[pointy, weighty ] =feglqdl (ngly) ;  %  quadrature  rule  for  y-axis 

quadrature  for  two-dimension 

for  intx=l:nglx  %  quadrature  in  x-axis 

point2  (intx, 1) =pointx ( intx) ; 
weight 2 (intx, 1 ) =weightx (intx)  ; 
end 

for  inty=l:ngly  %  quadrature  in  y-axis 

point2  (inty, 2) =pointy (inty) ; 
weight 2 ( inty , 2 ) =weighty (inty)  ; 
end 

-  end  of  feglqd2.m  - 
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o\°  o\o  o\o  o\o  o\o  o\o  o\°  o\o  o\o  o\°  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  o\o  hh  o\°  o\o 


Q.  glqd2 


feglqd2 .m 


unction 


[point 2 , weight 2 ] =glqd2 (nglx, ngly) 


Purpose : 

determine  the  integration  points  and  weighting  coefficients 
of  Gauss-Legendre  quadrature  for  two-dimensional  integration 

Synopsis : 

[point 2 , weight 2 ] =feglqd2 (nglx, ngly) 

Variable  Description: 

nglx  -  number  of  integration  points  in  the  x-axis 
ngly  -  number  of  integration  points  in  the  y-axis 
point2  -  vector  containing  integration  points 
weight2  -  vector  containing  weighting  coefficients 


determine  the  largest  one  between  nglx  and  ngly 

if  nglx  >  ngly 
ngl=nglx; 
else 

ngl=ngly ; 

end 

initialization 


point2=zeros (ngl , 2 ) ; 
weight2  =  zeros (ngl,  2)  ; 

find  corresponding  integration 

[pointx, weightx] =feglqdl (nglx)  ; 
[pointy, weighty] =feglqdl (ngly)  ; 

quadrature  for  two-dimension 

for  intx=l:nglx 

point2  (intx, 1) =pointx ( intx) ; 
weight 2 (intx, 1 ) =weightx (intx)  ; 
end 

for  inty=l:ngly 

point2  (inty, 2) =pointy (inty)  ; 
weight 2 ( inty , 2 ) =weighty (inty)  ; 
end 


points  and  weights 

%  quadrature  rule  for  x-axis 
%  quadrature  rule  for  y-axis 

%  quadrature  in  x-axis 

%  quadrature  in  y-axis 


end  of  glqd2.m 
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